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Section  I 


INTRODUCTION 


An  objective  of  the  present  research  program  was  to  develop  a  finite-element  based 
procedure  for  analysis  of  free-edge  de  la  mi  nation  specimens  using  through-the-thickness 
elements  and  including  Ixith  siretching  and  Itending  effects.  This  necessitates  the  use  of 
multilayer  plate  ihet tries  which  can  simultaneously  consider  bending  and  stretching. 
Multilayer  plate  theories  have  been  developed  using  assumptions  on  displacements  or  on 
strewjes.  I'he  former  class  of  theories  may  be  classified  into  two  groups,  viz., 

1.  Theories  based  upon  assumed  variation  of  displacement  as  a  polynomial  in  the 
transverse  coordinate  over  the  entire  thickness. 

2.  Theories  based  on  assumption  of  piecewise  linear  variation  of  displacements 
over  the  thickness  with  ‘nodes'  at  the  interfacial  surfaces. 

The  first  group  of  the  displacement  theories  has  been  found  to  be  inadequate  for  repre¬ 
sentation  of  behavior  of  composite  laminates  where  the  material  properties  along  fiber 
are  significantly  different  from  those  in  the  directions  across  fiber.  The  second  group 
of  theories,  called  the  discrete  laminate  theories,  apparently  is  a  better  candidate  for 
further  consideration.  These  discrete  laminate  theories  have  been  described  by  Srinivas 
[l]  and  Sun  and  Whitney  [2],  among  others,  and  solution  schemes  have  been  proposed. 
However,  for  arbitrary  geometry  and  a  large  number  of  layers  one  has  to  resort  to 
numerical  procedures.  The  finite  element  method  has  been  extensively  used  for  the 
analysis  of  plates  (  e.g.  Reddy  [3]  and  [4],  Davis  [5]  among  others). 

Davis  and  Mawenya  [5]  proposed  a  general  finite  element  formulation  using 
quadratic,  isoparametric,  multilayer  plate  elements  which  allowed  layers  to  deform 
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locally  with  no  restriction  imposed  on  the  relative  properties  ol'  the  constituent  layers. 
Since  this  formulation  considers  the  transverse  shear  deformation  in  all  layers,  it  is 
applicable  to  any  arbitrarily  layered  plate.  However,  the  stresses  are  discontinuous  across 
the  interfaces. 

As  a  starting  point  in  the  present  research  program,  Davis  and  Mawenya’s  approach 
was  used  to  develop  a  finite  element  solution  for  analysis  of  laminated  plates.  Davis 
and  Mawenya  did  not  give  details  of  the  theoretical  and  numerical  formulation  they 
used.  In  this  report,  the  theory  is  restated  in  variational  form,  and  an  extension  of 
the  general  variational  theory  is  specialized  for  implementation  in  a  finite  element 
computer  program.  Its  effectiveness  for  analysis  of  plates  and  also  its  inadequucv  in 
modelling  stresses  in  free-edge  delamination  specimens  are  noted. 

Section  11  contains  a  summary  of  the  equations  governing  bending  and  stretching 
of  laminated  plates.  In  this  section,  the  kinematic,  equilibrium,  and  constitutive 
equations  for  a  lamina  are  derived  based  on  an  extension  of  Mindlin’s  theory  of  plates. 
The  displacement  field  is  assumed  such  that  the  rotation  of  each  lamina  is  an 
additional  variable  independent  of  transverse  deflection.  The  set  of  coupled  field 
equations  and  interlaminar  continuity  conditions  is  written  as  a  self-adjoint  matrix  of 
operators.  Consistent  boundary  conditions  are  identified,  and  a  general  variational 
formulation  for  the  purpose  of  finite  element  approximation  to  the  problem  is 
developed.  Section  III  discusses  the  finite  element  formulation  and  computer 
implementation  of  the  theory.  Some  illustrative  examples  and  comparisons  of  results 
against  some  alternative  solution  schemes  are  discussed  in  section  IV.  This  section  also 
includes  application  of  the  procedure  developed  to  analyse  a  multi-layer  free-edge 
delamination  specimen. 
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Section  11 


BENDING  AND  STRETCHING  OF  LAMINATED  PLATES 

2.1  INTRODUCTION 

StructuraJ  elements  comptused  of  an  arbitrary  number  of  orlhotropic  layers  can  be 
approsi mated  by  finite  element  procedures.  In  these  composite  elements,  each  layer  may 
have  different  thickness  and/or  elastic  properties  and  different  orientation  of  axes  of 
material  symmetry.  In  the  following,  we  summarise  the  governing  equations  for 
bending  of  plates  based  on  the  following  assumptions. 

1.  Loads  are  carried  primarily  by  bending  and  stretching  of  the  plate. 

2.  Sliding  of  one  layer  past  another  is  impossible. 

3.  Plane  sections  normal  to  the  undeformed  surface  of  each  layer  remain  plane 
but  not  necessarily  normal  in  the  deformed  configuration  and  the  in-plane  dis¬ 
placements  vary  linearly  over  each  layer. 

4.  The  transverse  displacement  is  independent  of  the  transverse  coordinate  i.e.  the 
transverse  strain  vanishes. 

5.  Deformations  and  rotations  are  small  and  the  material  is  linear  elastic  so  that 
the  linear  theory  of  elasticity  is  applicable. 

Since  the  neutral  axis  is  a  priori  unknown,  bending  and  stretching  are  coupled  with 
respect  to  an  arbitrary  plane  of  reference. 

With  the  above  assumptions,  using  the  rectangular  Cartesian  reference  frame,  the 
kinematic  field  variables  consist  of  three  displacement  components  for  an  arbitrary 
point  on  a  reference  surface  defined  by  a  constant  value  of  the  transverse  coordinate 
in  the  reference  (undeformed)  configuration  along  with  the  values  of  the  rotations  of 
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segments  of  the  ray  along  the  transverse  axis  through  this  point  defined  by  intersection 
with  interlaminar  surfaces.  The  number  of  field  variables,  therefore,  is  2m+3  where 
m  is  the  total  number  of  layers. 

2.2  EQUATIONS  GOVERNING  BENDING  AND  STRETCHING  OF  LAMINATED 
PLATES 

2.2.1  Introduction 

The  generalized  equilibrium  equations  represent  the  integral,  over  the  thickness,  of 
the  three-dimensional  equilibrium  equations  and  of  the  first  moment  of  the  equilibrium 
equation.  The  constitutive  equations  are  .stated  for  a  linear  elastic  monoclinic  material. 
Tor  implementation  in  a  Ritz  type  finite  element  approximation  procedure,  the  problem 
is  formulated  as  a  set  of  self-adjoint  field  equations  with  consistent  boundary  operators. 
The  index  notation  is  used  throughout.  Latin  indices  take  on  the  range  of  values  1,  2, 
and  3  whereas  Greek  indices  take  values  1  and  2.  Subscripts  following  a  subscripted 
comma  denote  partial  differentiation  with  respect  to  the  coordinates  defined  by  the 
subscripts.  Summation  on  repeated  indices  is  implied  except  where  indicated  otherwise. 
A  pair  of  indices  within  parentheses  denotes  the  symmetric  part  of  the  tensor  described 
by  the  subscripts  and  a  single  super-  or  subscript  within  parentheses  denotes  'no  sum' 
on  that  index. 

The  actual  displacement  vector  at  any  point  is  a  function  of  the  coordinates  (x,) 
of  the  plate.  Assumption  of  transverse  displacement  being  independent  of  the 
transverse  coordinate  i,  makes  it  a  function  of  (x„)  only. 
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2.2.2  Kinematics 


{*om|X)nents  of  tl)e  displacement  vector  Tor  each  lamina,  assuming  linear  variation 
ol  inplane  displacement  over  the  thickness  ol  the  lamina  (I'ig.  1),  can  be  stated  in  the 
form: 

fi) 

)  )  =  uix  ) 

St  s  Of  S  €k 


Here  a  rectangular  (Cartesian  frame  of  reference  is  used,  uj**  are  components  of  the 


displacement  vector,  and  vj*’  are  the  'in plane'  displacements  at  the  kth  interface.  0!** 


are  the  components  of  the  rotation  of  the  kth  layer  in  the  a  — 3  planes.  For  infinites- 


imal  elastic  deformation  i.e.,  — ^-<<1,  the  strain-displacement  relationship 

0^, 


is: 


e  =  4-(ti  .-fu  )  =  li, 

Therefore, 

JO  _  U)  JtMO  _  JO  .  JO  (*> 
(*)  _  1  r-.U)  j  jWK 


(2) 


U1  ^  J  +0W1J 

oj  0  3,0  ^  o 


C3) 


where 


=  =iCv<'Uv‘'>) 

OjS  2 


(4) 


(O  _  ,(*)  _  1  fAi)  ,AiU 
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Figure  1:  Configuration  of  a  multilayer  composite  plate. 
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2.2.3  Equilibrium  Equations 

The  three-dimensional  equations  of  motion  including  body  forces  are; 

+  =  (5) 

where  cr,j  are  components  of  the  Cauchy  stress  tensor  and  /,  those  of  the  body  force 
per  unit  volume.  Equation  (5)  is  difficult  to  satisfy  exactly.  To  eliminate  the 
dependence  on  the  coordinate  Xj,  (5)  can  be  restated  in  the  form 


■* 

/,  'n  ,  j/'  u)  ‘inw  (o  VI  .  (o  -  „  , 

+  -p  u  Xx,  )  dXj  =0.  n  =  0,l,... 


(6) 


As  an  approximation,  in  general.  (6)  is  enforced  only  for  n=0,l.  Higher  order  theories 
\xould  use  higher  values  on  n  as  uell.  Equation(6)  for  n=0  is: 

a.  j-or  i=l,2,  setting  following  (l): 


P 


or 


where 


k 

.jii}  r  u  ■  j  (n 

/V  I  CT  .dx^ 

ikp  f  of*  3 

J'  /pdx'p 


(7) 


(8) 


(9) 


and 


f  p'‘'(Ux"')dx'‘' 


(10) 
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b.  For  1=3. 


or 


where 


r,  (o  ,  (/)  ,  (/)  "lo^  .  (o  ,, 

J  K3.n+<^33.^'*-A-P  ^^^3  =‘* 


’  +(o-^*  -  (T )+  M*’  -  F^'' ^  w'* W 

*'Qf,Of  33  33  3 


(11) 


r  (n  .  in 

G,.  =  J 


(12) 


,.0)  f  yJ-O  J  '-‘t 

^3  =  /.I 

*T> 

For  n=l,  i-e^  taking  the  first  moment,  considering  i=l,2. 

_(*)  .  _t*)  .  >(* )  -<*) :  P )  J* )/  J* )  ^  K  r  I  >1 J  J  - 


/ 


J.  '.y  '  ^'*  V  v“  ' -  M 

^'^«W+‘^<.3,3+4  "P  ~P  ^^3  ^^.^^^3  ^"^^3 


(1.3) 


or 


where 


»  +  li  " q'^  ’  -  v“ '  0“ ^=0 

I  a3  a  o 


_  r  j-t) 

-  J  4  ^3 


(14) 


(15) 


(16) 


and 


/*>=]'p“>u‘*V<> 


(17) 


In  the  approximate  theory  considered,  the  first  moment  equation  for  i=3  is  ignored. 


S 


To  eJitninate  the  time  derivatives  in  (7),  (11),  and  (14)  and  to  include  the  initial 
conditions,  convolution  of  both  sides  of  the  equations  with,  gCt)=t,  a  function  of  time 
is  performed  [8]  and  [9].  Hquations  (7),  (11),  and  (14)  can  then  be  written  in  the 


following  form: 

(18) 

(19) 

;a'^’=o 

(20) 

where 

(*)  represents  the  convolution  product,  livaluating  the  convolution 

integrals  for 

terms 

involving  time  derivatives 

g*  J  g(.i-r}w^\x^,T)dT 

*0 

=  f  (t-r) 

=  —  tw^*^(x^,0)  +  w'*^(x^2)  —  w^^\xp,0) 

(21) 

where 

and  are,  respectively,  the  initial  conditions  for  transverse 

velocity  and  displacement.  Similarly, 

(22) 

g*  ^  =  -  (  0“ \x^,0)  -  ^(i-^.O)  +  0^'  \x^4) 

(23) 

Let 

(24) 

r<*>  =  -p^^\  _  tv;;’(jp,0)-v‘;\3:^.0)J-  -  f 

(25) 

Z‘*>  =  - 1  0“\x^,O)-0|^*Vx^,O)] 

(26) 

Substituting  (24)  through  (26)  into  (18)  through  (20): 


9 


(27) 


(28) 


(29) 


Eqs.  (27)  through  (29)  are  the  spatial  equilibrium  equations  for  the  motion  of  each 
lamina  of  the  laminated  plate  in  terms  of  laminar  force  resultants. 

2.2.4  Constitutive  Equations 

2.2.4. 1  Stress-Strain  Relations  for  Linear  Elastic  Materials 

A  material  is  said  to  be  ideally  elastic  if  the  material  completely  recovers  its 
original  shape  upon  the  removal  of  the  forces  causing  the  deformation.  The  generalized 
Hooke's  law  relates  the  nine  components  of  stress  and  the  nine  components  of  strain 
by  a  linear  relation.  Assuming  an  initially  unstressed  reference  configuration,  and  a 
rectangular  Cartesian  reference  frame,  this  can  be  expressed  as: 

cr,,  =  £  e,.  (30) 

t/kf  kt 

where  are  components  of  fourth  rank  isothermal  elasticity  tensor.  Owing  to 

symmetry  of  tr,j,  i.e.,  in  the  absence  of  body  couples,  we  have: 

Furthermore  since  €,.^  =  €^^.1 

E...=E..,.  (32) 

tjk(  ijft 

If  a  strain  energy  function  exists  then 


(33) 


For  a  monoclinic  system  with  two  orthogonal  planes  of  symmetry  the  stress-strain 
relations  for  any  layer  k  in  the  reduced  form  with  respect  to  a  global  plane  of 
reference  (second  rank  tensors  written  in  vector  form)  are: 


1(1 


„(0 

*1111 

p(0 

*'1122 

^IJ33 

0 

0 

,.(0 

^'1112 

1 

^1  1 

CT  *  • 

li^l 

j/i} 

0 

0 

lO 

2i 

2211 

2222 

2233 

*'2212 

22 

J*’ 

*^33 

*'3311 

^0 

*'3322 

^0 

*^3333 

0 

0 

^3312 

(O 

'33 

^23 

0 

0 

0 

Uk) 

^2323 

p(*) 

^2313 

0 

■^*23 

^13 

0 

0 

0 

p(0 

*1323 

Ut) 

*1313 

0 

2€‘*^ 

'^®13 

°^I2 

p(0 

*1211 

jy<*) 

^1222 

p(*) 

*1233 

0 

0 

r<*) 

*'  1212 

^'12 

Hquation  (34)  can  be  written  in  the  indicial  form  as; 

''ofl  ^'nP33*3J 


U  )  _  -„(*)  JO 

cr  ,  =  2/j  ,  , 

o3  f»3y3  y3 


,J*)  Jt)  ,  ^l)  Jt) 


^33  ^33y6^y6 ^3333*33 


Solving  (37)  for  and  substituting  into  (35) 
U)  .  W*)  _«) 


where 


and 


U)  _  -jjfO  U)  ,  W*) 

tt/3  ^a0yt  €>6  '''  ^*^33  *^33 


pX*) 

Wi)  (*)  '^gg33  i7<*) 

aPyt  —  ^cifiyi  ~  J*)  33y6 

■^3333 


riU) 

Wi)  _  ■^0^33 
^0^33  _tO 

*3333 


Substituting  (3)  into  (38),  (36),  and  (37); 

(^)  _  t/*!  4.  j. 

afl  ^a$ySeyf>  ^afiybXs  *^o0  ■*■  ■“0,^33  ^33 

<r'“  =  £'*;  ,(w“’  +  *“') 

(Sf3  0f3y3  ,<sr 


cr  -  =  J^!^l  + 


,(0 

33 


e; ;  +  f;:.  ,  +  £**] 


33y6*3  '^yfi *3333  ®33 


(34) 


(35) 

(36) 

(37) 

(38) 


(39) 


(40) 


'33y6  y6 

We  note  here  that  the  assumption  of  «,*’  constant  over  the  thickness  of  each 


(41) 

(42) 

(43) 

lamina. 
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implies  u‘/’Cj:,)  =  w“'(x,,),  €'*'=<).  This  would  result  in  the  terms  containing  €,*’ 
dropping  from  (43).  We  also  note  that  this  would  mean  that  (37)  and  (35)  would 
not  contain  therefore,  there  would  be  no  question  of  eliminating  between 

these  two  equations.  However,  most  plate  theories  suffer  from  this  defect  viz.  e,*’  is 
eliminated  in  the  constitutive  equations  (35)  through  (37)  but  is  set  equal  to  zero  in 
the  equation  for  cr,,  (43). 


2.2.4.2  Constitutive  Relations  for  Force  and  Moment  Resultants 

Sulwtituting  from  (.35)  into  (i»)  and  (15)  iind  currying  out  the  integration,  the 
following  equations  are  obtained: 


^u)  _  hi  .  „(()  (*) 

<»/3  a/lyb^yh  a/lyli^yfi 


_  o<*)  f*)  j.  f*) 

~  ^«fiy6ey6'*'^„fiy6f<yi, 


where 


fM)  «)  (*)  ,  hjJt)  ,  , 


,(*) 

3 


(44) 

(45) 


(46) 


Substituting  from  (42)  into  (12),  and  carrying  out  the  integration  gives: 


^  fr 


,U)  /  (Jt)  ,  Aih 


&3y3  ,o 


where 


G‘*3  3  =  t 

v3y3  I  of3>3  3 


(47) 


(48) 


The  quantities  V<*’ ,  (2l*’ .  nIa  are  the  'laminar'  resultants.  It  is  to  be  noticed  that 
whereas  corresponding  quantities  for  the  entire  laminate  are  obtained  by  simple  addition 
of  d*’  and  is  the  total  moment  including  contribution  due  to  ,  i.e., 


e.  =  Zc?! 


(49) 


i-i 


0 


j 


12 


(50) 


* 


It 


I  J 


M 


(/£) 


(t-J 


(51) 


where  is  the  distance  of  the  center  of  the  kth  layer  from  the  point  of  application 
of  the  resultant 


2*2.5  The  Interlaminar  Continuity  Equations 

tor  lire  coniinuity  of  i factions  and  displacements  to  be  satisfied 


Hi ) 


a  .  =  cr  j 
1*,^  \ 


a  [ ' 


UJ  ,  ^  Af)  (JtM) 

V  +t,4>  =  V 

II  A  ■  (1  lA 


(52) 

(53) 

(54) 


iind 


a 


I ) 
33 


where 


(55) 


-  components  of  the  shearing  stress  at  the  top  of  the  kth  layer. 

“  components  of  the  shearing  stress  at  the  bottom  of  the  (k+l)th  layer. 
=  components  of  the  in-plane  displacement  for  the  kth  interface. 

<i)'A  =  components  of  the  rotation  of  the  kth  layer. 


2.2.6  Summary  of  Field  Equations  in  Convolution  Form 

For  self-adjointness  of  the  set  of  operators  consisting  of  those  appearing  in  the 
kinematics  equations,  the  constitutive  relations,  the  equilibrium  equations,  and  the 
continuity  equations,  noting  that  the  equilibrium  equations  had  to  be  transformed  (27) 
through  (29),  it  is  necessary  to  express  the  remaining  field  equations  in  convolution 
form  as  well.  i.e., 


1.3 


a.  I’quilibrium  1-quaticms  1(27)  through  (29)): 


r  f o-,,  -<^,,1  )  +  «*  -f' ,1  -  vv  +  =0 

o  *  0  5  o3  o  o  a 

b.  Kinematic  Relations  ((2)  through  (4)) 

s*e:i  =  jg*(v;;>vj;;) 

-fc  o.  U  )  j,  /  (/  )  I  j  )\ 

2g  €  1  =  g  +<^  J 

c.  ("onstituTive  liquations  ((44),  (45)  and  (47)) 

*  wU) _ *  M  _L 

”  -T  ^o^y6^y£ 

=  g*G^3^3C2€^3) 

d.  Continuity  Equations  ((52)  through  (55)) 

*  'Ki  )  ^  ^  ^  ■'f  ^  ^ 


*  ii}  *  (J<M) 

=  g*w 

*  +(x-)  *  — 

^0^33  =  ^^^^33 


The  field  operators,  for  layer  k,  in  self-adjoint  form  are: 


(56) 

(57) 

(58) 

(59) 

(60) 

161 ) 

(62) 

(63) 

(64) 

(65) 

(66) 

(67) 

(68) 
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% 


t 


0 

-a:''' 

0 

0 

0 

0 

0 

0 

g* 

0 

U) 

0 

0 

0 

0 

-8*^  t 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

g*I^ 

0 

0 

-g* 

[Af*’  = 

0 

0 

0 

it) 

-g*D^Py6  g* 

0 

0 

0 

0 

0 

0 

g* 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

_p(0 

0 

0 

0 

0 

0 

0 

0 

0 

O  ^QjyS 

g* 

0 

0 

0 

-g* 

0 

0 

g* 

®  aa 

0 

where 


L  =  1(8^  4-+S 

2  V"  9)3  yH  Qq 


The  domain  of  the  operators  is  the  direct-sum  space,  under  the  convolution  product,  of 
the  spaces  consisting  of  admissible  in  that 

order.  The  operators  in  (69)  are  the  ones  in  (57),  (62),  (59),  (58),  (63),  (60),  (56).  (64). 
and  (61)  excluding  those  associated  with  interlaminar  tractions.  The  displacement 
continuity  (66)  and  (67)  can  be  directly  incorporated  in  the  system  of  equations.  At 
the  same  time,  addition  of  operators  associated  with  the  interlaminar  traction  completes 
the  set  of  field  equations  for  the  system.  Traction  continuity  ((65)  and  (68))  is 
implicitly  satisfied  by  using  only  one  set  of  tractions  for  each  interface.  The  system 
set  of  field  equations  then  has  the  form: 


p''y  lo]  [d^) 

(0] 


{0} 

{Pf^ 

{0} 


(70) 


(0] 


[01  [d’^^ 

[c'-'  r 


lo-'"’-*’) 


{0} 
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or,  symbolically, 


[X]{y}  =  U} 


Here 


[5^'Y  = 


g*  0  0  0  0  0  0  0 

0000  00^00 


^*  000000  00 
0  00000  -^00 


i  (OiJ  ^  ui  ,(/)  .U)  ill  U>  ^  Hr  ,,(<) 


(71) 


(72) 

(7.1) 


\cr  ]  =  <or_^,  .cr,,  > 


33 


(F‘y  =  <-g*/-’'"-r'‘',o.o,-^-’;'-z'''.o,o.-g*y.;''-x'‘'.o.o 


,/n  „u) 

<1 


Jt)  .At  ' 


,.U  I 


(74) 

(75) 


{()}'  =  <{),()> 


(Vh) 


{o-'®y  =  <g^tr^3\o,0.0,0.0.^cr33’,0,0> 

=  <-rcrr3^lo,o.-r^g*<<;'\o.o.r<^^^ 


(77) 

(78) 


[0]  = 


0  0 
0  0 


(79) 


Continuity  of  tractions  is  ensured  explicitly  by  using  the  interfacial  traction  as  the 
field  variable  [cr'*T  in  the  manner  expressed  by  (70),  i.e.,  cr*^*’  =  and  cr'^d 

appear  as  field  variables.  The  operators  and  their  adjoints  [v9^*T,[C**r 

represent  coupling  between  field  equations  for  the  layers  and  the  continuity  of 
interlayer  displacements.  Explicitly,  for  the  interface  between  the  kth  and  the  Ck+lXh 
layer;  these  have  the  form 
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f 


f 


[Af^^ 


0  0  0 
0  0  0  0  0 


loj 


F 
0 
0 

0 

0 

0 

o 

0 

0  0  O  O  0 

0  jc*  0  0  0 
*:* 
o 
0 
0 

0 

0 

0 

D 

0 

0 

0 


0 

0 

0 

0 

0 

0 

g’ 

0 

0 

o 

0 

o 

0 

0 

0 

0 

u 

0 

0 

0 

0 


(01 


y 

tii) 


*0/3 


M 


rt*) 


2€’ 


ik) 


€’ 


-g*  0  O  0  0 
0  0  0  0  0 


0  0  0  0 
0  0  0 


^€Ei 

J*> 

^33 

'y 

t*is 

A*+l) 

Jk*l) 

V0 


I*4l) 

i* 

^*a3 


(*) 


»  u 


-g'C^**— 2**^ 


i> 


i) 

0 

■  K  0„  Z 


o 
o 

0 
0 
0 
0 


2.2,7  Self -ad  joint  ness  of  the  Operator  Matrix 

For  the  operator  matrix  [X]  to  be  self-adjoint,  a  sufficient  condition  is  that  the 
elements  of  X  satisfy  the  relationship 

<  /,  Boundary  terms 

Where  <  ,  >jt  is  a  bilinear  mapping  over  the  space  of  functions  defined  over  the 
region  R.  The  elements  of  X  satisfy  this  requirement  in  the  sense  of  inner  product 
defined  as 

</.«>=  f  /(£^8(x)dJi 

•5- 

Specifically,  if  the  operator  matrix  /4'*'  is  self-adjoint  and  ([^‘T, 

constitute  adjoint  pairs,  IX]  is  self-adjoint.  (Considering  /4,3  and  for  the  kih  layer. 
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based  on  Green's  theorem  iKreyszig  [1979]): 


<  >  =— g*vY^  > 


+  <  "L7  ’Ir  **'’“'  >  s*  ’1,  >  ,„, 


Similarly  for  and 


+  <  ^p  -8*  Kp  \  >  > 

■  t 


+  <  <p^p.r  >^u, 


iUK. 


to  */ 


and  for  ^4^  and  A^ji 


.»  yt) 


O  'O 


+ <  Q7%*g*('^y  >^*,+ < 


ii*' 


(80) 


(81) 


(82) 


The  remaining  elements  of  [A?*^  constitute  algebraic  operators  which  are  self-adjoint  or 
consist  of  adjoint  pairs.  The  adjointness  of  and  [//*’y,  and  of  [C**!  and  [C^^T  is 

obvious. 


2.2.8  Consistent  Boundary  Conditions 

Referring  to  Sandhu  [l975l  consistent  boundary  cx)nditions  for  the  problem  are: 


an 

on 

s'c'^q" 

an 

(8.1) 

on 

5*“ 

4 


» 
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% 


f 


f 


jCO  An 


W  w  =  g^v?*  T)^  on  S 


.in 

b 


where  C**’,C4**.C^*\C‘/’,C4*^,  and  C*/*  are  the  consistent  boundary?  operators  and 


I  Z  J  4  So 

•  £  J  4  >  u  ^ 


lixpticitlv, 

f'l  -f',!  =f^s  =-n 


(Corresponding  internal  jump  discontinuity  conditions  are: 
g*id;^N^^lr  =-g*igf)  on 

=-g*(«''3’)  on  5<3*/ 

^(^OgU)).  ^  -g.(g’f)  on 

=  g*(g'“'’)  on  s';/ 

g*{C^  <f>^)  =  g*ig  ^)  on 

U)y  */  'UK  rtU) 

g*CC(,  W  )  =  gHg  ^  )  on  S^. 


(84) 


Here  S',^\S*^\  ...  ,5^^’  represent  surfaces  imbedded  in  the  interior  of  the  region  R.  A 
prime  over  any  quantity  denotes  the  jump  in  that  quantity  across  the  surface  of 
interest,  e.g. 


/'  =  /"-/" 

where  ±  denote  the  two  sides  of  the  interior  surface.  Quantities  g’J*\  i-1,  2,  6, 

denote  the  specified  values  of  the  jump  discontinuities. 
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2.3  VARIATIONAL  FORMULATION  OF  THE  PROBLEM 


2.3.1  The  Basic  Variational  Formulation 

Based  on  (A .20),  the  self-adjoint  form  of  the  field  equations  given  in  (70)  the 
boundary  conditions  (83),  and  the  jump  discontinuity  conditions  (84),  the  governing 
function  for  the  variational  principle  is: 

N 


i-1 
A-J 


+  M  cr  >r,j+2l<y  ,c  cr  >  , 


(2.85) 


k^\ 

N-\ 


i‘2 


k-X 

—  2  <  >  —2  <  >  +  Boundary  terms 

Substituting  (70),  (83),  and  (84)  into  (85),  the  explicit  form  of  the  function  including 
the  boundary  and  the  discontinuity  conditions  is: 

N 


O  =  ]^{  <  -  y**  V  +  g*  p  -  r'*  VI*'  +  2g*  ^ > 

+  <  , -g*  ^ ^'11- 


,(*) 


(2.86) 


ot0yti  yh  j^*) 


+  <  V**  —  +e*e^*^> 


+  <  VI*' +  g*  -  g*  +  2g»  +  2Z<*^> 

+  <  »v‘*\  -  p'**  g*  +  2g*  F**' +  2x‘*'> 


Jf'*’ 


r^U). 


.U) 
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N~l 


•^L^<  \  >g*<^^3  >/*)+<  ^.  >y*)+<^  -^^<^33  >■/*)* 


Jt-1 


*-2 

A^-I 


/-I 


j^> 


-2<vl‘^g*<y'>^. -2<>v‘'\g*a3‘3‘’>^, 

+  Z{  <  ^l*^  +  2  i^l^2  V>  s(*> -^  <  ^ 

*-i  ‘  * 

+  <  w'*’ ,  fi  +  2  0,  V > +  <  Ql*^ ?*<  ci** w'*’ -  2  V>5tt) > 

+£{<  vl*\g^(df  y*j)'  +  2J^'<*^7ip>^,,+  <  <*’,«*C(C'2*’v‘*V-2«'!.*%p)>^*, 

t-i  *'  * 

+  <  .^'“ ,  |-(  (o'"  m")'  +  2  At*’  r,  )>,„+<  V“ ,  (C'*’  «“>)'  -  2  «■!*’  Tlj)  > 

Here  we  have  used  the  notation  for  o-q***’  as  explicitly  stated  in  (52)  and  (55). 

Theorem  The  Gateaux  differential  of  the  function  ft  defined  by  (86)  vanishes  if 
and  only  if  (70)  along  with  (83)  and  (84)  are  satisfied. 
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Proof. 

Let  =  be  an  admissible  state 

corresponding  to  the  set  of  field  variables 

/„(*)  _  j  J*)  J*)  A7<*)  o-‘*H 

Gateaux  differential  of  the  governing  function  (86)  along  the  path  B,  provided  the 
limit  exists,  is 

S.nCu)  =  +  +  (2-87) 


A-  I 


+  <  V  *  ’ ,  ^7t^l  „  -  F.  P*'  ^  -  f'*'  V'  ’  > 


+  <  «L" .  8*  *‘A’  -  8"  ’.j  V*  *“5  ^ 

+  <  35*" ,  -  g*  v“!„  >  ^,,  +  <  wl" ,  W.*i  -  8*7^^,  >  y,, 

+  <?“’.  g”  -  8*  el  -  VI*’  +  2  g-  g"’  +  2  Z"'  >  ^„ 

+  <  ■*■“’. g‘»l*6.,-g<*’-2i“’<‘’-/“’<*’  >^,, 

+  <  <r8*ML"-8*£‘'’.,v6'“-S'"“U,i''“  >,Ki 

''aP’8^  8^^  aPv6^>6  8^^ 

™  a/S  *  8^  '^oe  8^  *^(<».P)  ^  /« ■*'  ^  -^0.^  ’  ^  Jf(*) 

+  <  ^ d^l  -  +  2  g*  +  2  > y*)  +  < 

+  <  28".g>Q"-g“G<*’„3,3(2€")  >^„+  <  2€‘*>,g<*’-g>C'*’,3,,(2j"’)  >^„ 

+  <2'",g.(20-8-*l‘’-8'»“>,... 

+  <e">,g^2£'“)-g.?;;>-g.<'>  ,„) 
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N-l 


*-i 


jpt*) 


Jk) 


jfe-2 


+  £«<'’,-?•  (t;;;’  >  ^.,  +  <  w"’ ,  -  j*  cr-,‘,  > 


JV-J 


A-l 


+ <  <j’.s*<'’+<.8*<‘’-s^r‘’>.,.,+ <  <'’.w‘’-i*'<'‘"'>y.,) 


-  2  <  <‘  ’ .  J*  <^.■'3'  ’  >  ,,  -  2  <  »"’ .  S*  4'  ’  >  ,. 


-  2  <  l/,"' .  -  s*  ■r"’  >  -  2  <  ?■;■’ ,  -  (3,  »•  cr:;«  >  _,„  -  2  <  w*"' ,  -  *•  > 


AF 


+  £(<<*’.  -  w"  7,3  +  2  Xr“  71,  )>  3,,,  -  <  v" ,  J*  7?!%  > 

*-l  *  ' 

^7 


"  "t  's“> 

+  <€’•«*<-  O, + 2  *“>s  )>  ^,  -  < «“ .  S'  71,  > 

+  <  3<‘^«'<*l‘’7l3-2«“’7l,)  >,,.,+  <  «L*^s'»L'’71,  >,1., 

+  <  W“’.f'(-e“’71.  +  2(3“’71.)>^„-  <  3*“ 71.>^, 

+  <Sl‘’.S'(»“’71,-2ft“>7lJ>^„+  <Q“',jn5“’Tl.>^,) 

+  h  <  '“•«*<  -'<»>'  »3.«  -  <  >3... 

t-l  ■'  >' 


+  <  71, -2(<' 71,  )>  +  <  >3... 

*2# 
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+  <  7^<“ ,  WI'Vt),  -  2  («")■  7)j)  >  ^„  +  <  m“  ,  s*  (®“  ’)■  Tlj  > 

+  <  w“’ ,  -  (e‘*’)'  ii. + 2  ca.)'  i),)  >  jU)  ol*’)  »i.  >  jij) 

+  <S1‘\«*((«'‘V7,.-2C«<‘Vt).)>^„+  <Q'*’,S*(w'‘V7,,>^^.,) 

Similar  to  equations  (80)  through  (82).  the  following  relationships  hold 

<  >^.l  =  -  <  ^ 

<  ■*“.«*<«  >y..  =  -  <  <  ■*!.'’•**  O.  >sU. 

<  >,,.,  -  -  <  >^..+ >^., 

^  ^  y*)  =  —  <  ^  jj(*)  ■*■  ^  ^s**' 

<Q'‘’.f*«^’>y«  =  -  <e“’ .**''“’>^.,+  >j,., 

I 

where  5**^  =  Boundary  of  U  internal  surfaces  in  Using  these  to  eliminate 

2l*i  ^’'0"^  Gateaux  differential,  one  can  write 


JV 


8,  n(u)  =  2  £{  <  p‘*\ f*  -  V*’ ’  +  g* 


(2.88) 


*•1 


+  <  *1“ .  f  "1"  -  «*  *■ -  «■  ^“’.er»  "'i’ » 

s'.*’  •  *•  "Lm  -  «•  o'.  +  J*  c“  +  z“’  > 


+  < 


^  *^a0  ’  ^  ^  ^  ^  o$yb  ^  j/*) 

■*■  ^  Op’S* ''<1)3  ^^(0,0) 
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+  <  \  g*  Ql‘,i  -  ’  w''  ’  +  g*  4 


+2  <?'*\(,g*o-;;3*^>^,  + 


AT-l 


£-] 


/-z 

A^-1 


^2  D<<3Vvr+(,ff*c-^^c  >y*)+<^3;-s*>-  -^^  >/*> 


jt-i 


-2  <  <“.^*^^,3*  >3^-2  < 

+2  £{  <  9‘„*\  g^  -  y„*p  t)p  +  i^‘*p  )  >  3tt) 


Jt-l 


■^<^^\g<-  ^ >  5(*) 


+  <  R^*\  g*(  -  n„  +  0,  'n,)> 


si*’ 


^2  2:{  <  <*\g^  -cn‘:;)'  t,, + (O 

*-l  *' 
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+  <Jl‘“,j^(<‘’)'^,-(0")'V>s;;i 

+  <»‘*>,g^-(Q")'^,  +  ((2jT,.)>^„ 

+  <it‘',*^(w“’)'T,.-W“>r-iiJ>  ,., 

*^61 

Linearity  of  the  bilinear  mapping  <  .  >  implies  that  S^nCu)  vanishes  if  all  of 

equations  (70),  (83),  and  (84)  are  stilisfietl.  Also,  iiondegeneraieness  of  the  bilinear 

mapping  <  ,  >  implies  that  the  Gateaux  differential  in  equation  (88)  would  vanish, 

along  an  arbitrary  path  {«},  only  if  all  the  field  and  continuity  equations  and  tlie 

corresponding  boundary  conditions  are  satisfied. 


2.3.2  Extended  Variational  Formulation 

The  solution  of  (70)  must  belong  to  the  admissible  domain  of  the  operators  for  the 
functional  (86)  to  be  meaningful.  Using  (A.21)  along  with  (A. 22)  relaxation  of  the 
order  of  differentiability  of  either  component  of  the  following  pairs  is  possible,  viz., 
(Qa  or  wX  (J^ap  or  However,  relaxation  of  differentiability  could 

not  be  done  on  both  components  of  a  pair.  Some  possibilities  are  indicated  by  (80) 
throughd  (82).  Using  these  equations  to  eliminate  one  or  the  other  of  the  adjoint  pair 
of  operators  results  in  reduction  of  the  differentiability  requirement  on  some  of  the 
admissible  field  variables.  This  provides  a  basis  for  extension  of  the  variational  princi¬ 
ple  to  a  domain  where  the  differentiability  requirements  are  selectively  relaxed.  These 
extended  variational  formulations  also  provide  the  basis  for  certain  useful  specializations 
to  reduce  the  number  of  field  variables.  To  develop  these  extensions,  it  is  convenient 
to  rearrange  the  terms  in  (86)  to  write  the  governing  function  as: 


fl (u)  =  2  £{  <  ’ ,  f  Ff  +  r ''  S  , ,  +  <  0'^  \  g*  »  > 


I- I 
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N 


+  *^^"<>^.18  >=(*)+< ’'a  •  > 


^-J 


+  <  /^^  —  F*  > 

or^  *  ®  of^  ©  a^yfi  y6  &  »^yfe  >5 

^  ^  » 5  6  ^  ^ik} 


+  <0^  ,fM^^^^-g*Q^  -R  v„  -/  0^  > 


y*’ 


+  <  ff^* * , £*  Af** ^  —  e*  £** '  e‘' *  —  s*  *  K  * ^  > 


J.  ^  it  t')  _iJ.U)  ^ 


+  <  w‘*\ ?*  -  />'**  w‘* ’  > 


(n,fl)  '■^jp(*) 

jj(») 

+  <Ql*'.r(2€':i)-^0l*»-j^w^>>  } 


^f-i 


y*)' 
_+(*) . 


(89) 


+£«''".-*•  >  ^.. + <  '*'*’  ■-!-<  >  y.,  I 

*-2 

Af-1 

-2<v^  ,g*cr^3  >^,-2<>v  .y*cr33  >^ 

'*‘^‘*''^0  '8*  *^03  ^yAr)'^2<0|^  •*A'8*‘^o3  ^  yw)  ■*■  2  ^ ’8*  *^33  ^  J^N) 


+  £{  <  vl‘\  f*(  -  Af T)p  +  2  v‘*  ^  -  2  (»;;  * )  > 


^•1 


si*' 


+  <0'‘’,s*<  - 1,,  +  2  A“; T,J  >  ^,.,  +  <  m'‘1  T), ,g*< «'" - 2 «“>)  > 
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's  *^6 

+  h<  V'  +2fJ^l'i)'T)^)>^*,+  <  Op.g*(  Cv';>)'-2(<*’)') 

*-1  >'  ■’' 

+  <'#l‘’.«*<-("“iV'+^'*l«>'’'*  >si«+  ^  M'‘l-^,’f<<<i>'l'y-'^^o>T'>"»^» 

+  <  .g*(  -  (c;jS  j'  +  2  (Ol^’r  n„)  +  <  d\ .  g"(  (w‘‘>)'  -  2  V )  > ,,  ,i 


STliP  1 

For  exam  pie,  eiiminaiing  in  (89)  by  using  (8f)),  the  domain  of  N'*^  is  extended 

from  C*  to  (7,  and  the  function  has  the  form 

n ,  M  =  2  £i<  v",g*f'j’  +  y<‘>  >^.,+  <  *"’,j*g“’+z'‘’ >^., 

jt-1 

+  <w‘*’,g*f‘3*‘  +  X“*>^,) 

+  T  {  <  v'*\  -  v‘*’  -  je'*^  0'*’  > 

*  C  ^  or 


rXk} 


"*"  ^*0^'  «  "  a^yl^yb  o  "  o^yS^yS  ^  jf(*) 

^  ^  ^  "a^'5  ®oP  S  ^(or.a)  j,(*) 

+  <0^  ~R  -/  0„ 

OfflyfiS®  aPyS'^yS 


■*■  ^aP  ^  Aat./5)  ^  jp(0 


+  <2e2.rQ!f’-?*G"U3<2<P>,, 


.(*)  /„ 


,(*> 


(9(0 
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A'-l 


+  L’<^  -«  '^aS  .S*(^33  >  ^*)  > 


i*-! 

+  L<  <  -  g*  o-;;* >  >  +<„,(*>,_  g*  o--*  > } 


S-l 


+ z«  ‘^r>*’.s*v'‘’+<,**«L‘’-**vr'’  >^.,+ < 


,«)  *  U+l) 


,<l)  *  (1)^ 

-2<^  ,1*0-^,  >j,,-2<w  ,g*ff33  >^, 


+  2  <  .  g^  CT;«f  >  >  ,,  +  2  <  0-’ , g>  a;'f  >  ,,  +  2  <  w'-'’’ . g-  > 


AN) 


C.V) 


J-VI 


'^2 


+  •  2  «*  >  i,)  +  <  7J^ ,  2g*(  V 

+  <^1*\ r(  -  +  2  V  > /o  +  <  <^1*  -  2  <*»)  > 


+  <  >v'*>.g^  -Clf'T,^  +  20^11,)  +  <  Q1*'t,^ .«*(  2 

N 


+  D  <  v",2^/>'“,,  >  +  <  <^,,2^  V'»-<‘>)  >^., 

k~i  *»' 


+  <«“’ ,  j<<  -  Af  ■"  T,, + 2  *•;;>  7),)  >  ^.,  +  <  T,, ,  jn  *■!*’  -  2  <'  ’’  > 

+  <  »“’,s*(  -Q'l‘\,  +  2  2’"\)  >^„  +  <  Q“ti,  2  «'“’)>  J„) 


29 


STEP  2 


Eliminating  using  (81 )  leads  to 

*‘l 


/-I 


_l  ^  -»  xl*)  _*  -X 


+  2<N^^^  £*e^^— > 


4  JO 


■*■  ^  ’  S*  ttflyS  ^yh  ^  ^  ft^yS  ^  j^^k') 


+  2  <  *  S^  'f„p  *^(o.i9}  ^  j^t> 


,(*) 


,{*) 


JV-l 


(91) 


*-» 


I  ^  U)  *  “(jfe)  _1  ^  ^  'W  I 


*-2 

^-1 


jU)  *  a-Hi) 


Hi)  *_U) 


*-i 


-2<vl‘\g-crJ3‘'>  -2<>v'".g*c7;;’>^, 
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+ 2  <  + 2  <  € '  ■ '  .V  + 2  < 


N 


+  £(  <  v‘‘’.2**J»l>,  >^.,+  <  Nl>,,2r(v"-«“’)  > 

i*  I  * 

+  <0".2**a1>,>^.,+  <«L>,.2r<0L*’-O> 


4*^ 


4« 


+  <  w“> ,  -  q"  „.  +  2  6"  11„)  >  ,.,  +  <  qI*’  71. ,  j*(  w'*'  -  2  «'*’ )  >  ,„  ) 


-lU). 


5'*> 


+  Zl  <  v“’,2j*A’«7,,  >  +  <  A/"t,j.2s*(  v“’-c’“>)  > 

j{  J  '  H 


4f> 


+  <  -Q'“’71.  +  2!3'"t|.)  >^„+  <  S“’7l.,*^»'''‘’-2»'“’)  ) 


STEP  3 

Eliminating  Q[,*^  in  (91)  using  (82)  results  in: 


*-i 


+  <iv“\j*27‘3’  +  2f“’>^„) 


0  *  o 


i-1 


+  <«***  — o*J9***  > 

*e*  «  ^  oS>6*y6  S  ®  ttjSye'^ye 


+  2<W^*^  F*  —  ff*  > 


!>(*)„(*)  rWl.Ul. 

+  <0^  ,-je  -/  0^ 


+  <  >,c*> 

+  2  < 
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(92) 


JH\ 


A^-l 


.  ^  f  ^  U)  *  +U)  ^  ^  -  _*— j  \ 


.u)  ^  *U). 


i-1 


+  I«  v".-*-o-r,*’  >^.,+  <  w“'.-|-cr;*  >^„) 


,U) 


it-a 

N-l 


+  L^‘^°’«3  +*/g  “S  ''c  +  -s  > 


■  U)  .  (/•!) 


(')  ..  (O  i.  ('•!) 


<-1 


-2<v^  >^,-2<w  .g^o-,,  >^, 

_L  ^  (A’)  _*  _+<A)  ^  .  .(A')  ,  .  *{m  .  .  T  ^  „,(.v)  .  -  „'(.V)  . 

'?**^a3  ^  jf(A‘) ^  ^ ’'*'^°^33  ^  jp'A'i 


AT 


1-1 


+  Z*  <  ''"•2s*i»‘>,  >^„+  <  A/"t,,.2j^v"'-0“')  > 


U)  a(*>' 


s(*) '  - '  ''p-‘fi '  «  »  '  4*' 


+  <*“.2«*itr<.*’Tij  >^„+  <  Af'^Tls, 2*^ *“’-«“’)  >^„ 


+  <  »’"’.2j»fi''’T).  >^,1+  <  Q“’i1„.2*^»’“’-*'‘’)  >jU,  > 


JV 


+  Z(  <  v‘*\  2  g*  />■"  7,j  >  ,„  +  <  w'-;  T), ,  2  j*<  v'; '  -  o'- '  )■  > 


(i£) 


U)  M)^ 


4*’ 


„(*) 


+  <  >v“\2g*(2'l*S„  >^>+  <  Ql%..2g*(  >,*, } 


4*> 
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2.3.3  Some  Specialization 

The  function  defined  by  (92)  has  no  differentiability  constraints  on  the  stress 
resultants  function  also  leads  to  certain  useful  specializations  by 

identically  satisfying  some  of  the  field  equations.  Assuming  that  in  (92),  continuity  of 
is  identically  satisfied  i.e.  normal  displacement  is  restricted  to  be  constant  through 
the  thickness,  i.e.,  tv**’  =  w  for  all  k  (67), 


fii 


/•  \ 


(93) 
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+  2  <  .  ?*  >  ,.^•.  +  2  <  <f>T. ' J.  5*  +  2  <  .  +  g*  cr;';'  > 


+  Li  <  v“^  2  r  np  >  +  <  T)^ .  2 v‘;>  -  ’ )  > 


i-l  * 

jU)  «  * 


A 


+ <>v“\2rors„  >^.,+<(2::^„.2^(>v<^>-^"')>^,,,  i 

+  £{  <  vl-\2g^;^'l>,>  +  <  A;l>^,.2g^.r-0r)'  > 

/-I  "  ■’' 

+  .2g*( > 

*1/  '^4^ 


If  the  kinematic  relations  (59)  through  (61)  are  identically  satisfied,  equation  (93)  leads 


to: 


n ,  (u)  =  2  £{  <  vj,*’ .  g*  +  <  0'*^ .  g*  +  d^'>> 

*-l 

+  <>v'*\g*  7^3*4 

N 


,tjt> 


jt-i 


^  ^  a^y6  ^y6  ®  ^  a^y6  ^  ^ fc) 


+  <  0^*  \  -  ie‘*^  -  /* 

^Or  Or  0/ 

,  ^  U)  *  dU)  <^) 

+  ^  Jf  —  H  .  ^  . 


nt*) 


^  D'^*2  *  y^U)  U)  ^ 

^  ^  (^^y&^yS  6  ^  or^y6  ^y£  ^  /fr) 


+  <  2  el'] .  - ^ g"»„3^3  (2  el'])  > ^*,  +  <  w . -  A''' 4>  } 


C94) 
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i-1 
N 

+  T,^<  '''t’ <  "■  .-r  "j3-'j,„i 


+  L«  "'I'j' V  + 1,  s*  «!*’  -g*  v“'  ■’  +  >  +<<>■*>.*•  w“>-  j*  w'""  >  ) 


+  2 <  v“’ ,g* <'3'”  > ,,„  +  2  <  j'<'f  > _^,,,  +  2 <»“>,  +  j* cr;'f  > 

^  £(  <  v«,2g.;«l>,  >  +  <  <;„,.2gK  /'-«“')> 

*-l  ' 

+  <^“' .  2g»  ,,,  +  <  v‘; ,,, .  2*^  «“>-<*’)  >^„ 

+  <  w“’,2g*S“’Tl„  >^„+  <  Q“’r,„,2j<(  ».“>-«'*>)  >^„  ) 

+  D  <  vr,2g.O,  >  +  <  > 

*-l  *' 

+  <«“’.2J*A*>„^  >  +  <  Jh“;t)  2j«<0“-$«>)-  > 

+  <  *<'>)'  } 

If  the  kinematic  boundary  conditions  in  (83)  and  (84)  are  identically  satisfied,  equation 
(94)  reduces  to 

n,W  -  2£i<  v“',j*y«+y'‘'  >^„+  <<‘’.s‘c“+2'‘’>^., 


+  y  { <  V**  \  -  p'*  ^  -  p'*  '  > 

t-l 


jit) 
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Let 


i_  ^  U)  *  Ai)  (O  _*  D«)  U»  ^ 

'*‘‘^^0^’  aflvfc^vE  rvflvfi'^vE  ^ 


yh 


yit 


+  <0^  ,-R  v^  -/  0^ 


*  "  ofivE^vS  6  ^  «e>8"v6  ^ 


+  <  2e‘^*^,-g*G'*'„3^3(2<;)  >  ^*,+  <w.-P^‘''w> 


oPyE'^yS 
-.(*) 


.U)^ 

alyJ  '">3^ 


(95) 


A'“i 


X*l 

A' 


.U) 


*{n 


tn  .  'UK 


+  Li  < -r .- g*  f* ‘-a'.  >  ,u., ) 


i-2 
A- 1 


.  tr*  I  ^  +^i)  *  (i)  ,  *  *  jUJ  *  UMJ.  *  ^  *.J*^  **-^  X 


i(0  j. .!**<) 


<«)  j,..«)  *..U*0. 


i-l 


-  2  <  vl*\  g*  a:;‘S -  2  <  ^ 


fi 


»3  gl 

jfm 

(*)  **(*)„  ^  .  ^J.U)  ^  a(*)„  ^  ,  ^_.M) 


+  2<V^  ,g*ff„3  >  ^Af) ^ ^or  '*JV^*^o3  ^  j^JV) ^ ^  *^^*’^33  ^jf<.v> 


+2 D  <  C-r  J^;>^  >^*)+  ^ >s(*)+  <  ^  .5*61^.  >s.*) » 

*-j  *  ’  ^ 

+  2  £{  <  v“\g*  >^^,+  <0';\g»  J^’l'^T)^  >^„+  <  w<‘\2g»0t^.  >,UI  ^ 


*-i 


^3-0^33  °^33 


Then,  if  the  constitutive  relations  (62)  through  (64)  are  identically  satisfied,  Cljiu) 
specializes  to 


n,(u)  =  22;{<v';\g»F‘j’  +  y<*>>  ,+  <0l'\g*G^"’  +  z'*S 


,(*) 


*-i 


+  <w.ff‘‘'  +  x'‘'>^„) 
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N 


+  ZI-  <  <S'S* 

Jt-i 


N-l 


+  Z  ^  J  .  ?*  >  y*)  +  <  '  • '/  >  y*)  +  <>^ .  g*  0-?f  >  ) 


x-1 

.V 


+  Z<  <  ''o  ^  •  ■  g*  ^03  ’  >  ,r*.  +  <  W .-  g*  0-/3  >  } 


i-2 

A  I 


-  2  <  v‘ '  \  g*  4* » >  ^,  +  2  <  vl^  \  g*  <r 

+  2 <  g*  +  2 <  w.g* P3 


N 


+2ZU  <  w.g*o<.ii,>^^*)} 


+  2£{  <  v“  \  g*  T,^  >  5(«  +  <  <^1*^  •  «*  >  s<*)  +  <  «  .  f*  <3'„  } 

*-l  »'  ^S( 


ir  the  interlaminar  displacement  continuity.  (66),  is  identically  satisfied,  (96)  reduces  to: 


8  =  2  £{  <  v^*\ g*  +  y***  >,*,+  <  g*  +  z'*’  > 


*-i 


y*) 


+  <  K-g*  ’  -  /e‘* '  v*'  ’  -  ‘  '  > 
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(97) 


-2<  >^.  +2<  >^„,, 


JV 


+  2£{  <  >  ,*,+  >^,*,+  <  «'.r!5,T7„  > 

,.,  Si  s,  ij 

+2£i  >,j*>+  <’^'•8*0',. T)„  >  } 

*'l-  *’■!*  ■’t. 


Here,  vj,*’  is  no  longer  independent  for  k=  1,  2 .  n.  Purilierniore,  if  siirfuce  slieiir 

tractions,  surface  shear  couples,  body  forces,  and  b<xly  couples  are  neglected.  i.e.. 
=  =  =  =  0  then  (97)  reduces  to: 

n,(u)  =  2  2;(<  v“’y*>  >  ,„+  <  >_^.,+  <  w,y‘>  >^„) 


4-1 

N 


+ D  -  <  < .  **  C)  >  ,<«  -  <  "“i  ■**■**„  > 


*-i 


J?*** 


j.  ^  ^  n(*  1  d(*1  r'*'  > 

+  <0^  ,~g*Q^  ~R  ^  -/  >^u) 


/*> 


,(*)J 


(98) 


+  2  <w,g*i>, 


+ >^*,+ <0f  ,g*  >^*>} 

4-1  *  3  S 


;v 


+  22^{  <  >  +  <  W,g^(3'^T,^  >  } 


U)  »  A.U), 


t-i 
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Explicitly,  replacing,  using  the  continuity  of  the  displacements  across  the  interface, 
by  v”’  and  0“’,  U,  is  written  as: 


N 


i-i 


"  ,o<“)  =  2  <  ’■I”  •  +  <  Z  >,.«+<  >,«)+<»’.  z‘*’  >  ) 


N 


i-1 


-  Z<  <  <^  «* "!  ”» >  y.)  +  <  •  **  Z > 


i-i 


i-1 


y*) 


+  <  .  g*  <^>(^,6)  >  +  <  G,  . g*  w „  > 


+  2  <  >^ui  +  2  <  >y*, 


j^\nf  '^Q  '  V  ^pr  jptJt'J 

lU)  ^Ai)' 


y*) 


+  <  2  <  >y*) 


k-l 

..(I)  o(4)  .  ,(/) 


HI 


Jl-1 


+  <Z‘,*”>'’Z^'^^ 


Hi 


i-1 


yo 


+  2  <  w.^P  > 
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/"J 


*  i-  i  1 


i-1 


4*^ 


(99) 


+  2£{<  >yo+  <  Z*itf  >y*) 

i-J  HI 

+  <  >^„+  <  w,g*0>.>^„ ) 

If  inertia  terms  are  ignored  i,e.  for  the  static  problem,  it  is  not  necessary  to  take 
convolution.  Then,  the  governing  function  is  the  functional: 

n  ,  j(u)  =~Z^'^  ^np'S*^iy,h)  -^oP’g*  Z‘l‘^(y.«  ^  jpl*l 


i- 1 


I'.- 1 
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(1(X») 


+  2  <w.g*/>3  >j, 

+  2£«  vl‘>,**i»«r,,  >^.,+  <  > 

'  f-1 

+ 2  £{  <  v;." .  g*  T)^  >  ^,,. ,  +  <  £  ^ ,  ,  s* 

t  I 

+  ^  ^  ^  ^  ^  ,1  Ui  ^  ^^-Ul  ^ 


1^ 


If  the  physical  problem  does  not  have  line  loads  or  couples  applied  to  surfaces  in  the 
interior^  ^nd  vanish.  This  gives  the  governing  functional  as 


N 


n  ,  j  C«)  -  -  £{<  .g*  >  Jt)  +  ^(y^)  ^ 


k~\ 


+  <  >/..+  < a'"-e‘<^^  >^„+  <<'.s*Q'">^.,  > 

+  2<w.g*P3>j,  (101) 

This  specialization,  dropping  the  convolution  with  g(t),  was  used  to  set  up  the  Tinite 
element  approximation  discussed  in  the  next  section. 


ffk) 

0)  -«to. 
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Section  III 


FINITE  ELEMENT  FORMULATION 

3.1  INTRODUCTION 

The  Cinite  element  method  subdivides  a  given  region  K  in  an  n-dimensional 

euclidean  space  into  a  number  of  disjoint  open  subregion  (elements)  =  ] , 2 .... ,m 
such  that 

k  =  lim  IJ  W  (102) 

Here  a  superscripted  bar  over  a  quantity  denotes  its  'closure'  i.e. 

^  =  R"  [J  (103) 

where  Q/?  is  the  boundary  of  R*.  Disjointness  of  these  elements  implies 

R'’r|R^  =  0  if  eP^f  (104) 

A  set  of  nodal  points  in  O  defines  the  geometry  of  the  elements.  In  the  following 

section  a  discretization  of  the  domain  by  the  finite  element  method  is  presented.  The 

formulation  of  the  finite  element  is  based  on  the  variational  principle  governed  by  the 
specialized  functional  n,j. 

3.2  FINITE  ELEMENT  DISCRETIZATION 

Let  the  field  variables  at  any  point  within  an  element  be  represented  by: 

»(<'>) 

=  ))(<”)  (105) 

*v(:ep  =  UV} 
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Here  are  sets  of  interpolation  functions  relating  the  values  of  the 

corresponding  variable  at  the  nodal  points  to  an  arbitrary  point  within  the  elentent  and 
{V"’},  and  {W}  are  vectors  of  the  a  priori  unknown  values  of  .  w  at 

the  nodal  points.  Substitution  of  (105)  into  the  governing  functional  (101)  yields  the 
spatially  discreti7.ed  functional  as: 

m 

f*  I 


where,  far  N  layers  in  the  laminate: 

A 


'  I  i 

+  2  r {*;>  )dR'' 

Jf 

'i'  '■*  '■* 

+  2^*  die 


k-i 


2  f 

■  f  {/>///(<!>!;’}  d/e' 

I  O  ^  V  7 

[DH] 


f  ,{H  f ) 


iW} 


(106) 


dST  } 


+  2  {//..} /’,d/e‘’ 


Here  the  symbol  13  associated  with  H  denotes  appropriate  differentiations;  [yA'*T.  L^*’] 
and  tZ)‘*1  represent  constitutive  relations;  and 


1 


t 


f 
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A-;;’  =  f 

Pa  5 

=  f  0HjlA^‘){DJ//dR‘^ 

=  /  \DH^}{A^'^\DH/dR'- 
fT 

=  /  {/)//,]  [a" ’]{/>// 

I 

P' 

=  /*  lo//^)[/y'*]jz>y//rfA‘- 

/ui 

Ja) 

^!iic  =  / 

3») 

R  =  f  i//  )P,dR‘ 

h/,  j  V  i 


where 


iDHf  =0  // 
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y*y 


Ky 


0 


[Z)///  =0  Hi 


4Ky 


„r  r 
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[dhJ  = 


[V  = 


ll\  0 


0 


Using  the  above  definitions,  equation  (l()6)  can  be  written  as: 

A  J.  I 

n  = y  f  {4>'n 

c  <'  t  ) 


t  I 


i  1 


y'J 


#^1 


i-t 

+  {4>[;  {W}  + 

+  2  {H^}^  R 

vp 

Gateaux  differential  of  the  governing  functional  (107),  denoting  by 
path  of  variations  in  (Vj  ,{<!>„},  {w}  respectively,  is: 


N 


i-l 


S.ft  = 

u  e 


k~i 


-  2  £1  ir<J  >)’■  ( v';>) + >)’■  [at'^  ( £  r,  («;» 

1  i^h 

(v' '  ’) + z  wii"  iCi  z <*;i 


j-1 


1-1 


j-l 


+ ’  y  if'*’} + i^,;  ’}''  ’1 


(107) 
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it- J 


/  1 


+ lOd  <,  (<»  +  (£■, 

i  i  t  I 

+  {w} 


(108) 


+  2{R^^/i 

k/i 

1'he  vanishing  of  the  Gateaux  differential  S^fl  (KJS)  for  arbitrary  {I^J,  {$J, 
yields  the  following  equations. 


r- 1 


il“l 


i-l 


=  0 


(109) 


m  N 


r«  I  i:  -  J 


i-1 


(no) 


Jt-l  ?-j 


f-l 


1-1 


and 


Z<  -  lO  ('*'1  -  )+<'<..,)  I 


t-i 


=  0 


(in) 


Where  the  ^  denotes  the  "direct  stiffness”  addition  of  contributions  of  all  elements. 
#^1 

Equations  (109)  through  (111)  can  be  written  collectively  in  the  following  matrix 
from, 

t^]{t/}  =  {R) 

Here 


i/f]-  £iA'r 


c- 1 


45 


U}  = 


and 


lul  = 


where  n  is  the  total  number  of  nodal  points  in  tlie  system  and  we  identify: 

2, N*2 

3, JV+2 

4, N-*2 


[Kr= 


^.3 

., ..  /f,. 

^LA-.I 

K 

^22 

^23 

^24  •■  •■ 

^2j+1 

■•  ^2.  A' 

^3,A'4  I 
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^33 

K,,  ..  „ 

"  "  ^3,N 

1/ 

^3.JV4i 

K 

K,,  ..  „ 

"  ■■  ^4,N 

^4, a;.. 

K 

• 
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•  ; 

t+l^*  1 

*  * 
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^  NJ^+l 

K 

K  K 


^  AN2*A'4  2 


where 
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=  ,LiO 

^l 

+  Z'O 

/-iV-J 

f-N 
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=  B<'J 


h  J 
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=‘20  +  '2‘/-lZlO 

M 


/'  V  I 


^2.A-,.  =‘2l'^!i”l  +  '2'.V  .ZI<’J 


/  A' 


‘^4JVN2 
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Af 
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/-JV 
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and 


^A-.2.»<2=lO  +  fO 


{vY,  =  <  v'  ,...,0'  ...0"  > 


O  J 


UK  = 


up 


where  the  subscript  j  denotes  the  jth  node. 


(112) 


3.3  SOLUTION  PROCESS 

Since  the  stiffness  matrix  is  symmetric  and  banded,  only  the  above-diagonal  terms 
in  the  global  stiffness  matrix  are  stored  in  the  form  of  a  rectangular  matrix.  The 
first  column  of  this  matrix  is  the  diagonal  of  the  original  matrix.  The  load  vector  is 
set  up  according  to  (112).  After  the  global  stiffness  matrix  and  load  vectors  are 
assembled,  the  system  of  simultaneous  equations  is  solved  for  the  displacement  degrees 
of  freedom  using  the  Gaussian  elimination  process.  The  solution  for  nodal 
displacements  (S')  of  an  element  along  with  [5]  matrix  which  is  the  element 
St  rain -displacement  relation  are  used  to  calculate  the  stresses  at  the  desired  points  of 
the  element.  i.Cn 

n 

e  =  as'  =  7  5  8 

I  I  I 

/-I 

where  n  represents  number  of  nodal  points  in  the  element. 

(T  =  ^58' 
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and  contains  the  elastic  constants  of  the  specific  lamina. 


in  (106)  the  specific  [/i] 


is  denoted  by  [DHY 


Section  IV 


VERIFICATION  AND  APPLICATIONS 

4.1  PROBLEM  DESCRIPTION 

The  theory  discussed  in  Sections  11  and  III  was  implemented  in  a  computer 
program.  The  computer  program  was  written  in  I’ortran  and  designed  to  run  on  the 
!I1M3<)S1  maiivrranie. 

highl-tnxled  serendipity,  nine-noded  l.agrangian,  and  four-noded  Lagrangian  quadri¬ 
lateral  isoparametric  plate  elements  were  used  (Fig.  2).  The  nodal  degrees  of  freedom 
(the  number  of  unknowns  associated  with  each  nodal  point)  depend  upon  the  number 
of  layers  in  the  laminate.  If  the  number  of  layers  is  m,  then  the  degree  of  freedom 
is  (2m+3)  for  each  nodal  point.  This  consists  of  three  displacement  components  (u,v,w) 
for  a  reference  point  on  the  'through  the  thickness'  nodal  point  and  rotations  (0^.  0^) 
of  each  layer  to  describe  completely  the  deformed  geometry.  Two  versions  of  the  pro¬ 
gram  were  prepared.  The  first  was  an  "incore”  program  wherein  all  the  matrices 
were  generated  and  stored  in  core  and  the  solution  process  did  not  require  auxiliary 
storage.  This  version  was  adequate  for  solution  of  example  problems  and  verification 
of  the  genera]  approach.  Later,  a  version  was  developed  where  the  algebraic  equations 
were  assembled  and  stored  in  blocks.  TTiis  version  was  used  for  solution  of  a 
22-laminae  problem  using  the  CRAY-XMP  computer  at  the  Ohio  Supercomputer  Outer. 
This  implementation  was  verified  by  application  to  several  example  problems.  These 
examples  included; 

1.  A  simply  supported  square  sandwich  plate  made  of  isotropic/orthotropic  surface 

layers.  The  core  was  assumed  to  have  a  finite  shear  modulus  but  zero  elastic 
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modulus  in  extension. 


2.  1  Tee-edge  del  ami  nation  specimens 

a.  Angle-pl,v  l±45l 

b.  Cross-ply  [0/90\ 

The  results  were  compared  with  available  solutions.  The  effect  upon  the  accuracy  of 
the  results  of  mesh  refinement  and  .sublayers  subdivision  was  studied.  The  code  was 
used  to  determine  stresses  and  deformations  in  a  multi-ply  free-edge  delamination  speci¬ 
men. 


4.1.1  Analysis  of  Plates 

’I'iie  multilayer  plate  finite  element  method  was  used  to  ana)y7.e  a  three  layer, 
square,  simply  supported  sandwich  plate  uniformly  loaded  in  the  transverse  direction. 
The  geometrical  and  material  properties  of  the  plates  were  (the  top  and  bottom  layers 
are  denoted  by  subscripts  1  and  3,  the  core  by  2).  This  example  was  the  same  as 
used  by  [5] 

Plate  dimensions; 

Length  of  each  side  »  10  inches. 

Tydckness  Ij  and  of  surface  layers  =  0.028  in 

Thickness  of  core  1,  =  0.75  m 

Material  properties 

a.  Stiff  layers  are  isotropic  elastic. 

E,-E^~]0’  Vb/in 
^  =  >^3=  0-3 
Gj  =  3x10''  Vb/vn 

b.  Stiff  layers  are  orthotropic  elastic. 

E  ,  =  1 0^  Ib/in  ' 

1  I  ii 
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\ 


/i  =  A'  = 

y  I  y-^ 

a  =(;  =  1.875*10*'  Ib/in- 

A  vj  t,v3 

p  .  —  p  ,  =  0.3 

A  1  A  3 

C^,2  = 

=  1.2*10'*  Ib/in^ 

lauding 

g=l  Ib/Ui^ 

I7ue  to  double  svmmetry,  only  one  i)Uiirler  «r  the  plyte  was  analyzed  in  each  case 
Cl'ig.  3).  Wlien  eight-  and  nine-noded  elements  were  used  four  mesh  discretiza  tors 
viz.,  Clxl),  (2X2),  (4.\4),  and  (8X8)  were  used.  In  the  case  of  four-noded  element,  the 
discretization  for  the  quarter  of  the  plate  was  extended  to  a  finer  (16x16)  mesh. 
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A 


v‘''  =  0f  =  O 


(One  quarter  of  the  plate) 


Figure  3:  Mesh  configuration  and  boundary  conditions  for  a  quarter  plate. 


4.1.2  Analysis  of  Four-Ply  Free-Edge  Delamination  Specimens 

The  free-edge  delamination  specimen  was  treated  as  a  special  case  of  a  plate 
rig.(4).  The  boundary  conditions  considered  were:  fixed-fixed  (longitudinal  displacement 
specified)  at  the  longitudinal  ends  and  free-free  at  the  transverse  edges.  The  results 
for  four-ply  symmetric  laminate  with  stacking  [±45],  and  [0/90],  were  compared  with 
Pagano's  solution  [lO]  which  is  based  on  a  generalization  of  Reissner's  theory.  The 
comparison  covered  only  ,  and  because  Pagano's  solutions  were  available  only 

for  these  quantities. 

The  dimensions  were  specified  as  a'b  ^  1  and  b  =  8h,  where  a,  b,  and  h  are 
length,  width,  and  thickness  of  each  lamina.  res|wctively.  I'or  proper  comparison,  the 
material  properties  assumed  were  the  same  as  in  [lO]. 

=  20x10*’  psi 

=  2.1x10*  psi 

=  C,3  =  Gj,  =  .85x10*  psi 

V  ^  u  ~  ~  Oi21 

U  23 

Following  Pagano,  N»2  indicates  each  lamina  is  treated  as  a  single  layer,  and  N=6 

indicates  that  each  lamina  of  thickness  h  is  modeled  by  three  sublayers  with 

thicknesses  of  h/3.  For  N=)0  each  lamina  with  thickness  h  is  subdivided  into  five 
sublayers.  In  going  from  three  sublayers  to  five,  only  one  of  the  sublayers  was 

subdivided  into  three  new  sublayers  with  thickness  of  h/9. 

Two  different  types  of  -  mesh  refinement  were  considered.  In  both  types, 
refinement  was  carried  out  in  a  manner  such  that  the  nodal  points  of  the  previous, 

coarser  mesh  were  a  subset  of  the  finer  mesh.  Refinement  associated  with  longitudinal 
direction  was  performed  by  dividing  the  domain  in  the  longitudinal  direction  into  the 
stated  number  of  equal  length  elements.  In  the  case  of  the  transverse  direction,  the 


eo 


X  =  COrtsiarti  plane 


]''igure  4:  Configuration  for  laminated  coupon  under  uniaxial  tension. 
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refinement  was  in  the  form  shown  in  Figs.(5)  and  (6}i  to  ntinimize  the  ret^uired 
storage.  In  one  sequence  of  refinements,  shown  in  Fig.CS),  the  width  of  the  specimen 
was  first  divided  into  five  equal-width  elements.  Only  the  edge  elements  were  refined 
further,  each  time  the  edge  elements  being  subdivided  into  three  elements.  This  process 
was  continued.  In  the  other  sequence  illustrated  by  Fig.(6),  an  edge  element  was  used 
such  that  the  (y/b)  ratio  for  its  center  was  0.995.  The  remaining  interior  domain  was 
discretized  into  three  strips.  Refinement  over  the  thickness  was  done  on  the  central 
sublayer  while  changing  from  N=6  to  N=10,  Fig.(7c)  in  order  to  get  good  estimates  for 
T,,,  and  cr,,.  On  the  other  hand,  to  evaluate  t,.  at  the  interface,  discretization  was 
performed  on  the  sublayer  next  to  the  interface,  I'ig.CTd).  Because  of  symmetric  stacking 
sequence,  only  two  layers  were  used  in  the  analysis. 

4.1.3  Results  of  the  Analysis 
a.  Rectangular  Plates 

Figs.(8)  through  (ll)  show  a  comparison  of  the  approximate  solution  with  that  of 
the  sandwich  plate  theory  as  presented  in  ISj.  Fig.(8)  shows  the  convergence  of  the 
central  deflection  of  a  sandw'ich  plate  w'ith  isotropic  layers.  Results  using  the  Q8  and 
the  Q9  elements  converge  to  the  exact  solution  more  rapidly  than  those  from  the  O'* 
element.  TTie  results  obtained  for  one  element,  in  the  case  of  09  element,  or  four 
elements  in  the  case  of  Q8  elements,  are  superior  to  results  obuined  from  256  Cl6xl6) 

Q4  elements.  Fig .(9)  shows  that,  for  a  given  error  limit,  the  amount  of  time  needed 

/ 

for  Q8  and  09  elements  is  considerably  smaller  than  for  04  elements.  Results  for  a 
sandwich  plate  with  orthotropic  layers  are  presented  in  Figs.(lO)  and  Cl  l).  The  same 
observation,  can  be  made  as  in  the  case  of  isotropic  plates.  The  numerical  results  for 
isotropic  and  orthotropic  cases  are  presented  in  1’ables  (1)  and  (2),  respectively. 

The  program  was  also  used  to  solve  the  isotropic  problem  where  the  shear 
correction  factor  was  ignored.  The  maximum  central  deflection  for  an  8x8  mesh,  using 
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y 


I 


(b)  Nine  subdivision  across  the  width 


Figure  6:  Sequence  of  refinement  across  the  width  with  a  thinner  edge  ele¬ 

ment. 
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(b)  Six  sublayers  (N=6) 


Figure  7:  Sequence  df  refinements  over  the  thick niw;. 
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X  Q9  ELEMENT 
+  08  ELEMENT 
A  04  ELEMENT 


g  O  SERIES  SOLUTION 


Figure  8:  Number  of  subdomains  along  each  side  vs  central  deflection  of  a 

square  isotropic  simply  supported  plate  with  No.  of  elements  - 

N\ 
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X  Q9  ELEMENT 
+  Q8  ELEMENT 
A  04  ELEMENT 
O  SERIES  SOLUTION 


Figure  10\  Central  deflection  of  a  square  orthotropic  simply  supported 
sandwich  plate* 
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T<AU  J:  Numerical  values  for  maximum  deflection  and  the  cofrespondidg 

CPU  time  for  different  types  of  elements.(Isotropic  case) 


Number  of 
Elements 

Element  Type 

CPU  (Sec) 

Maximum  deflection 

1x1 

0* 

0.046 

0.000064402 

2x2 

CK 

0.154 

0.00036986 

4x4 

Q4 

0.605 

0.00059685 

8x8 

Q4 

2.918 

0.00069959 

16x16 

Q4 

19.897 

0.00072964 

1x1 

Q8 

0.232 

0X1006189 

2x2 

Q8 

0.934 

0.00073851 

4x4 

08 

4.492 

0.00073997 

8x8 

08 

29.377 

0.00074005 

1x1 

09 

0.287 

0X10073539 

2x2 

09 

1.221 

0.00074087 

4x4 

09 

6.516 

0.00074008 

8x8 

09 

50.446 

0.00074006 

Series  Solution 

0.00074 
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Tabu  2:  Numerical  values  for  maximum  deflection  and  the  corresponding 

OHJ  time  for  different  types  of  elements^Orthotropic  case) 


Number  of 
Elements 

Element  Type 

CPU  (Sec) 

Maximum  deOecdon 

1x1 

Q4 

0.046 

0.00026536 

2x2 

Q4 

0.146 

0.00066988 

4x4 

Q4 

0.S91 

0.0010218 

8x8 

Q4 

2.851 

0.0011671 

16x16 

Q4 

19.071 

0.0012078 

1x1 

Q8 

0.228 

0.0010439 

2x2 

Q8 

0.954 

0.0012194 

4x4 

Q8 

4.497 

0.0012215 

8x8 

08 

29.467 

0.0012216 

1x1 

09 

0.286 

0.0012140 

2x2 

09 

1.224 

0J0O12224 

4x4 

09 

6.545 

0.0012216 

8x8 

09 

50.267 

0.0012216 

Series  Solution 

0.00123 
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Q8  element^  was  =  OlKK)436()l.  'fhe  series  solution  for  same  problem  gave  almost 

identical  value*  =  0.(XK)43605.  I'herefore,  maximum  deflection  due  to  shear  for 

this  plate  would  be: 

W„  -  0.000740  -  0*0004360  =  0.000304 

^max 

b,  Four^Ply  Frce-Edge  Delamination  Specimen 

Values  of  0^3,  presented  in  the  subsequent  sections  were  determined  based  on  the 

constitutive  equations  (36)  and  (37).  Determination  of  these  stress  components  is  also 

possible  through  the  use  of  the  equilibrium  equations.  However*  inplanc  strain  is 
linear  over  un  element  and  discontinuous  at  the  nodal  points,  resulting  in  a 

discontinuous  set  of  in  plane  stresses.  Furthermore,  the  numerical  evaluation  of 
needs  one  numerical  differentiation,  and  needs  two  numerical  differentiation. 

Therefore*  the  use  of  numerical  differentiation  would  not  result  in  a  satisfactory 
estimate  unless  the  number  of  elements  in  the  y-direction  is  increased  to  a  point  where 
the  discontinuity  of  is  reduced  considerably.  This  would  result  in  a  very 

expensive  computational  analysis.  Therefore*  it  was  more  convenient  to  determine  <r^ 
from  the  constitutive  relations, 

I.  Angle  —  Laminate  [±  45X 

Since  the  stresses  were  determined  at  the  center  of  the  elements,  to  determine  the 
predicted  results  at  the  center  of  the  longitudinal  dimension  (x  =  L/2),  an  odd  number 
of  elements  were  used  when  using  Q8  or  09  element.  Because  thrac  elements  have 


midside  nodes*  displacements  were  also  directly  available  at  x  =  — • 

2 


In  the  case  of  04 


element*  two  analyses  had  to  be  carried  out*  To  obtain  the  solution  for  str^ses,  an 
odd  number  of  elements  were  used.  However*  as  for  the  QS  and  09  elements,  to 
obtain  the  displacement  results  at  (L/2)*  it  was  necessary  to  discretize  the  longitudinal 
direction  into  an  even  number  of  elements*  (odd  number  of  nodal  points)  so  that  a 
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set  or  nodal  points  would  be  located  at 
04  element 

Fig.  (12),  shows  plots  of  values  of  cr^  at  mid'surface  of  the  top  lamina  for  mesh 

refinement  in  the  y-direction.  The  number  of  elements  in  the  x-direction  was  kept 

constant  at  11  and  over  the  thickness,  each  layer  had  constant  rotation  about  the  x- 
and  y-axes.  As  observed  from  the  figure,  accuracy  of  the  axial  stress  is  unaffect^  by 
refinement  in  the  y  direction.  However,  the  values  do  improve  somewhat  with 
refinement  in  the  longitudinal  (x)  direction  (Figs.  13  and  14).  Fig.  (13)  shows  the 
effect  of  refinement  along  x  into  21  and  31  elements  with  the  number  of  elements 
along  y  kept  constant  at  17.  Fig.  Cl4)  shows  the  effect  of  refinement  along  x  for  9 
elements  in  the  y-direciion.  Fig.  (15)  shows  the  effect  of  thickness  refinement  on 
The  results  improve  significantly  with  refinement  near  the  free-edge.  The  same 
observation  can  be  made  for  in  Figs.  (16)  through  (18)  respectively.  The  calculated 
values  of  improve  slightly  with  refinement  in  the  x-direction  but  are  unaffected  by 
the  refinement  along  the  y  and  z-directions.  Accuracy  of  T„  is  not  affected  by 

refinement  in  the  longitudinal  or  the  transverse  directions,  Figs,  (19)  through  (21). 

Accuracy  of  the  stresses  at  the  free-edge  improves  significantly  with  refinement 

over  the  thickness,  i.e.  as  N  is  increased  from  2  to  6.  This  refinement  does  not 
introduce  additional  elements  but  increases  the  degree  of  freedom  at  each  of  the  nodal 
points  by  assuming  roution  to  be  constant  over  a  smaller  portion  of  the  thickness. 
Comparison  of  results  for  different  values  of  N  indicates  that  N-6,  in  comparison  to 
N=2  and  N-10,  gives  more  accurate  prediction  for  cr^  Fig,  (15).  Further  refinement  to 
N=10,  in  comparison  to  N«=2  and  N«6,  leads  to  better  estimation  of  Fig.  (22)  and 
T,.  Fig.  (23).  It  should  be  mentioned  that  thickness  refinement  does  not  influence  the 
accuracy  of  stresses  near  V/b=0.  Stresses  at  the  free-edge,  however,  are  greatly 
influenced. 


69 


70 


71 


72 


+04  ELEMENT  N=10  21X9 

A  04  ELEMENT  N=6  21X9 

0  04  ELEMENT  N=2  21X9 
o  -  pfiGRNO  N=6 

(\j 

• 

cn 

i 

o 

o 

1 

i 

1 

i 

I 

i 

i 

i 

i 

1 

i 

1 

1 

i 

i 

1 

1 

. 

r . 1 

i 

( 

» 

{ 

. i 

1 

i 

i 

j 

■■  • . 1 

i 

t 

~9" . 

$  1 

£  08  ’Z  09  • 

z 

jsd'i  01/ (tij  •A 

1 

1 

\ 

1 

E 

1 

+ 

i 

t 

O 

1 

1 

I 

! 

..  i 

^  -  j 

1 

; 

t 

i 

t 

[ 

1 

1 

1 

1 

..J 

3 

i 

5 

1  1 

1 

1  . .  , 

^  i 

;  t 

! 

i 

’ . o . . 

i 

. — 

j 

.  ...' 

i 

i 

[ 

\  I 

. 

j  1 

i 

E . 

. . . . 

. . . 

1 

! 

1 

i 

1 

r 

1 

i 

1  ,1 ,,  ^ 

1 

A  c 

b 

o 

1 

! 

\  H 

1 

\  1 

1 

t 

s 

V 

CsT 

o 

ca 

i 

I""" . . 

s 

. . 

i 

i 

i 

1 

1 

1 

E 

3 

Figur 

'.00  o'.  20  o'.  40  o'.  60  o'.  80  1 

Y/B 
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Figure  16:  XY-stress  at  midsurface  of  top  layer  with  refinement  in 

y-direction;  Angle-ply  specimen. 
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X\'-stre$s  at  midsurface  of  top  layer  with  Ttfinement  in 
x-direction;  Angle-ply  specimen. 
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ligure  20:  XZ-stress  at  z=h  with  refinement  in  x-dircclion:  Angle- ply  spec¬ 

imen. 
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Figure  21:  XZ-stress  at  z»h  with  refinement  in  x-direction;  Angle-ply  spec¬ 

imen. 
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Comparison  of  accuracy  for  dif  ferenl  types  of  elements 


Accuracy  of  the  stresses  calculated  using  different  types  of  elements  was  compared 
and  the  execution  time  for  each  case  noted.  Four,  eight,  and  nine-noded  elements  were 
used.  Two  different  mesh  sizes  viz.  9x9  and  21x9,  with  edge  elements  were  used 
with  four-noded  elements  (Fig.  6).  For  eight  and  nine-noded  elements  only  a  9x9 
mesh  with  'edge  elements',  i.e.  elements  along  the  edge  having  very  small  dimension  in 
the  y-direction,  was  used.  The  CPU  times  on  the  Cray  X-MP/28  for  this  set  of 
problems  are  shown  in  Table  (3).  Fig.  (24),  for  N=2,  shows  that  the  Q4  element  with 
the  21x9  mesh  predicts  the  axial  stress  at  the  free  edge  with  11.477o  error  as  compared 
to  the  08  element  with  ll.CWf-  and  Q?  element  with  10.99%.  For  N-6,  Fig.  (25),  the 
error  for  the  (J4  element  reduced  to  2.68%,  for  QS  to  2.22%,  and  for  09  to  2.12?  in 
Fig.  (26),  for  N-10,  the  error  for  0^  's  5.72?e  and  for  08  is  5.11%.  Throughout,  tne 
04  element  gives  significant  error  near  Y/b-0.  Comparison  of  inplane  shear  stress 
shows  that  the  results  at  the  free  edge  are  predicted  with  approximately  same  accuracy 
(Fig,  (27)  through  (29))  by  all  the  elements.  In  the  case  of  (Fig.  (30)  through 
(32)),  04  with  9x9  and  21x9  mesh  refinement  gave  practically  the  same  results 
throughout  the  width  of  the  specimen  showing  that  the  system  is  rather  insensitive  to 
refinement  along  the  length  of  the  specimen. 

Further,  09  gave  the  same  result  at  the  free  edge  but  were  more  accurate 

than  04.  In  Fig.  (33),  the  longitudinal  displacement  predicted  at  X-L/2  and  the  top 
surface  of  the  specimen  using  08  element  is  shown  for  N»2,  N-6,  and  N=10.  In  Fig. 
(34)  and  (35)  corresponding  results  are  shown  for  N-2  and  N-6  in  the  case  of  09  and 
04  elements.  It  is  seen,  Fig.  (33),  that  for  N-2  the  displacements  are  slightly 
overpredicted  at  the  free  edge  and  for  N=6  and  N-IO  displacements  are  slightly 
underpredicted.  Though  the  results  improve  with  refinement  from  N=2  to  N*=6,  further 
improvement  with  refinement  to  N=10  is  not  realized. 
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Table  3:  Comparison  of  CPU  time  on  X-M1V28  for  Q8,  and  Q9  ele¬ 

ments. 


N 

Element 

Mesh 

CPU  on  X-MP/28 

(SEC) 

9X9 

1.82 

o 

21X9 

4.479 

08 

9X9 

16.555 

09 

9X9 

28.181 

04 

9X9 

16.961 

a 

21X9 

40.279 

08 

9X9 

190.222 

09 

9X9 

323.116 

04 

9X9 

60.755 

10 

21X9 

145.421 

08 

9X9 

704.333 
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X  Q9  ELEMENT  N=2  9X9 
+  Q8  ELEMENT  N=2  9X9 
A  04  ELEMENT  N=2  21X9 
O  Q4  ELEMENT  N=2  9X9 
'  PflGflNO  N=6 


Figure  24:  X-stress  for  N=2,  Q4,  Q8,  and  09  elemenis. 
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Figure  25:  X-stress  for  N=6,  Q4,  08,  and  09  elements. 
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Figure  27:  XY -stress  for  N-2,  Q4^  08,  and  09  elements. 
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X  Q9  ELEMENT  N=6  9X9 
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A  Q4  ELEMENT  N=6  21X9 
O  04  ELEMENT  N=6  9X9 
-  PflGflNO  N=6 


Figure  28:  XY-stress  for  N=6,  Q4,  QS.  and  Q9  elements. 
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Figure  29:  XY-stress  for  N»10,  Q4,  Q8,  and  Q9  elements. 
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Figure  30:  XZ-stress  at  z=h  for  N=2,  0^.  Q8,  and  09  elements. 
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Figure  31:  XZ-stress  at  z»h  for  N=6,  Q4,  Q8,  and  Q9  elements. 
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Figure  32:  XZ-stress  at  z=h  for  N=10,  Q4,  Q8,  and  09  elements. 
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Figure  33i  Longitudinal  displacement  at  U(L/2,y,2h)  of  Q8  element,  N»2,6 
and  10. 
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A  Q9  ELEMENT  N=6  9X9 
O  Q9  ELEMENT  U=2  9X9 


^  "  PflGRNO  N=6 


Figure  34:  Longitudinal  displaceinent  at  LJ(L/2.y,2h)  of  Q9  element,  N=2, 

and  6. 
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Figure  35:  Longitudinal  displacement  at  U(L/2,y,2h)  of  Q4  element,  N=2, 

and  6. 
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U.  Cross' ply  laminaie  [O  '90], 

Values  of  t,..  and  O’.,  were  calculated  for  four-layer  cross-ply  [0/90],  using 
equations  (2.33)  and  (2.34)  as  was  done  in  the  case  of  t,.  for  the  angle-ply  laminate 
[±45], 

As  shown  in  Fig.  (36)  mesh  refinement  in  the  y-direction  does  not  lead  to  any 
noticeable  improvement  in  the  calculated  values  of  Fig.  (37)  indicates  that 

approximation  of  is  not  dependent  on  axial  refinement  either.  Results  obtained  for 
refinement  through  the  thickness,  Fig.  (38),  for  11x17  mesh,  17  being  the  number  of 
elements  along  the  y-direction,  indicate  that  this  refinement  profoundly  influences  the 
results.  However,  the  traction-free  boundary  condition  at  the  free  edge  is  not  satisfied 
and,  therefore,  the  solution  does  not  match  Pagano's  results.  The  calculated  valuer.  ,r 
CTy.  are  in  error  up  to  309o  even  with  refinement,  for  y/b  equal  to  0.8.  Near  the 
free-edge,  the  error  is  quite  large.  The  values  of  Cjj  determined  from  equation  (37) 
were  either  exactly  ?ero  or  close  to  zero.  This  is  clearly  wrong  and  represents  a  serious 
limitation  of  the  theory. 

4.1.4  Analysis  of  22-Layer  Free-Edge  Delamination  Specimen 

The  procedure  developed  was  applied  to  a  22-layer  coupon  with  fiber  orientation 
of  [( 25.5 /- 25.5  )s / 90],  which  was  previously  solved  by  Chang  [ll]  and  Dandan  [12]. 
The  laminate  width  and  ply  thickness  were  taken  as  1.0  inch  and  0.00505  inch, 
respectively.  The  material  properties  of  a  lamina  were  the  same  as  used  in  previous 
investigations  (ll],  [12]. 

=  19.26M0*  (psO 

=  £33  =  1.32*10^  ipsi) 

G,,  =  G,3  =  G33  =  0.83 »  10*'  (psi) 

^12  ~  *'1.1  “  ^2.1  ~ 


+  Q4  ELEMENT  N=2  11X17 
A  Q4  ELEMENT  N=2  11X13 
O  Q4  ELEMENT  N=2  11X9 


°  ■  PflGflNO  N=6 


Figure  36:  YZ-stress  at  z=h  with  mesh  refinement  in  y-direction;  Cross-pi y 

specimen. 
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Figure  37-.  YZ-stress  at  z-h  with  mesh  refinement  in  x-direciion;  Cross-ply 
specimen. 
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Figure  38i  YZ-stress  at  z*h  for  (11x17)  mesh  in  case  of  N^2,  N-6,  and 
N^IO;  Cross-ply  specimen. 
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IDue  to  the  symmetry  ol'  t!ie  laminate,  it  was  only  necessary  to  consider  11  lamina  in 
the  analysis.  Bjised  on  the  analysis  performed  on  the  four-layer  coupon  specimen,  it 

was  concluded  that  refinement  along  the  x-direction  improves  the  accuracy  of  the 
results  more  effectively  than  the  refinement  along  the  y-direction.  For  this  reason,  the 
22-Iayer  specimen  was  discretized  into  22  elements  in  the  x-direction  and  13  in  the 
y-direction.  All  elements  have  the  same  dimensions  in  the  x-direction.  However,  in 
the  y-direction,  the  edge  elements  have  a  width  of  t).(K)5  inch,  with  the  remaining  por¬ 
tion  discretized  shown  in  i'tg.  (6).  'I'he  formulations  in  [ll]  and  [12]  were  based  on  a 
two-dimensional  model  dependent  on  y-  and  /  ciKirdinates  only.  In  [12),  22  elements 
were  used  in  the  z-direction  and  14  in  the  y-direction.  In  the  present  investigation, 
the  discreiization  is  along  the  x-  and  y-coordi nates.  Therefore,  it  was  not  feasib;'"  to 
match  the  mesh  with  that  used  in  [ll]  and  [12]. 

To  compare  the  results  with  those  given  in  [12]  and  [11],  two  different  analyses 
were  performed.  In  [12]  the  results  are  given  at  the  center  of  each  layers.  In  the 

present  investigation,  components  are  calculated  at  the  centers  of  the  laminae,  and 

based  upon  the  effectiveness  of  the  refinement  over  the  thickness  noticed  in  the 
analysis  of  the  four-layer  specimen,  each  layer  is  divided  into  three  sublayers,  forming 
a  total  of  33  layers,  is  calculated  at  the  interface  in  the  present  investigation. 

However,  in  order  that  these  re.sulis  be  comparable  to  those  obtained  in  [12]  and  [11] 

at  the  center  of  the  lamina,  each  layer  must  be  subdivided  into  two  sublayers,  forming 
a  total  of  22  layers,  producing  the  necessary  interfaces.  An  average  axial  strain  of 
(0.95414*10"*)  was  applied,  which  is  the  same  as  applied  by  Chang  [ill 

Fig.  (39)  shows  the  cross-section  of  the  upper  half  of  the  symmetric  laminate  and 
defines  the  location  for  plotting  of  results.  Fig.  (40)  and  Fig.  (41)  show  the 
distribution  of  tr^  along  the  center  of  the  11th  and  the  5th  layer  for  different  (y/B) 
ratios.  The  results  obtained  from  FSDT  match  those  obtained  from  the  higher  order 


element  [Hl  and  the  axisymmetric  model  [12].  Similarly,  Fig.  (42)  shows  that  at 

the  center  of  the  11th  layer  can  be  obtained  accurately,  though  the  concentration  of 

the  stress  at  the  free-edge  does  not  match  that  obtained  by  the  higher  order  element 
[ll]  Distribution  of  tr,,  along  as  shown  in  Fig.  (43),  closely  follows  the  result 

given  in  [ll].  However,  the  stresses  determined  from  the  FSDT  are  in  considerable 
error  near  the  free  edge.  Figs.  (44)  to  (47)  show  the  plots  for  cr^  along  Rl,  R5,  R6, 
Rll,  and  Figs.  (48)  to  (50)  indicate  the  cr„  along  R5,  R6,  and  Rll.  It  is  evident  that 
both  er,.  and  cr^.  are  predicted  reasonably  well  by  the  FSDT  in  the  regions  y/5<<).7. 
However,  for  y'S>0.9,  the  predicted  stres.ses  for  a,,  even  result  in  signs  different 

from  those  in  [ll].  Results  for  er^.  do  not  satisfy  the  free-edge  stress  at  y/B  «  1.0. 
F’igures  (51)  through  (53)  show  the  through-the-thickness  distribution  of  ff..,  cr^.,  and 
<r,,  along  the  free-edge.  The  results  from  the  FSDT  do  not  agree  with  those  from  [llj 
and  [12]  and  are,  apparently,  quite  wrong.  Furthermore,  Figs.  (54)  to  (57)  show  the 
distribution  of  o'„  along  Rl,  R5,  R6,  and  Rll.  Tlie  stress  distribution  obtained  from 
the  present  approach  is  quite  different  from  that  given  by  [12]. 
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Figure  39:  Cross-section  of  the  upper  half  of  the  22-layer  coupon. 
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figure  40:  Y-siress  at  R1  for  22-layer  coupon. 
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Figure  42:  X-stress  at  R1  for  22-Jayer  coupon. 
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Figure  43:  XY-stress  at  Rll  for  22-layer  coupon. 
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Figure  46:  ZX-stress  at  k6  for  22-layer  coupon. 
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Figure  49:  YZ-stress  at  R6  for  22-la.ver  coupon. 
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Fifftre  SO".  YZ-stress  at  Rll  for  22-layer  coupon. 
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Figure  51:  Through  the  thickness  distribution  of  Z-stress  at  y/b-0.995  for  the 

22-layer  coupon. 
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Figure  52:  Through  the  thicknew  distribution  of  XZ-stress  at  y/b-0.995  for  th 

22-layer  coupon. 
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Figure  53:  Through  the  thickness  distribution  of  YZ-stress  at  y/b-0.995  for  the 

22-laycr  coupon. 
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Figure  54:  Z-stress  at  R1  for  22-Jayer  coupon. 
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Figure  55:  Z-stress  at  R5  for  22-Iayer  coupon. 
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Figure  56:  Z-stress  at  R6  for  22-layer  coupon. 


119 


22  LAYERS  ZZ-STRESS  RT  Rll 


^  ORNDRN 

§  O  23X13  N=22 

m 


Figure  57:  Z-stress  at  Rll  for  22-Jayer  coupon. 
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4.1^  Summary  and  Conclusions 

1.  The  08  and  09  elements  are  superior  to  the  04  element  for  the  analysis  of 
plates.  However,  for  analysis  of  coupons  under  uniform  extension,  the  per¬ 
formance  of  the  08  and  0?  elements  was  comparable  to  the  04  element. 
The  CPU  times  needed  for  the  0®  and  09  were  much  greater  than  required 
for  the  analysis  using  04  elements.  Thus,  the  Q4  element  appears  to  give  the 
better  combination  of  accuracy  and  economy  of  computational  effort  for  cou¬ 
pons. 

2.  To  improve  the  predicted  results  for  o-^  and  at  y/b=0,  refinement  along 
the  longitudinal  direction  seems  to  be  most  effective.  However,  the  solution 
at  the  free-edge  depended  upon  thickness  refinement.  The  best  results  at  the 
free-edge  for  (7^  were  for  N=6,  while  for  they  were  for  .NslO. 

3.  The  results  for  t„  and  did  not  improve  significantly  with  refinement  in 

X-  or  y-directions.  The  thickness  refinement  improved  the  predicted  results  at 
the  free-edge  for  but  did  not  satisfy  the  traction-free  boundary  condition 
at  the  free-edge  for  t„. 

4.  The  predicted  results  for  cr,  are  all  either  zero  or  close  to  zero  for  the  four- 
layer  delamination  specimens.  However,  for  the  22-layer  delamination  speci¬ 
men,  they  were  oscillatory  along  the  z-coordinate  and  quite  different  from  the 
results  obtained  by  Chang  [ll]  and  Dandan  [12]. 

5.  The  predicted  displacements  for  free-edge  delamination  specimens  are  close  to 
Pagano's  results.  For  N»2  the  displacement  results  are  slightly  overpredicted, 
and  for  N=6  and  N*10,  the  displacements  are  slightly  underpredicted. 
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Section  V 


DISCUSSION 

t 


In  the  present  research  program,  the  theoretical  studies  carried  out  have  included 
two  distinct  approaches  to  the  problem  of  stress  analysis  of  composite  laminates.  One 
approach  consisted  of  a  specialization  of  the  three-dimensional  elastostatics  theory  to  the 
ca.se  of  free-edge  deJaminution  specimens  under  uniform  axial  strain  in  which  the  stress 
field  is  independent  of  the  longitudinal  coordinate.  The  other  approach  consisted  of 
development  and  application  of  theories  of  laminated  plates  to  the  problem  of  fre;-edge 
delamination.  The  present  report  covers  study  of  the  applicability  of  the  existing 
"discreu  laminate  theory”  as  one  of  the  first  steps  in  the  research  program.  Thts 
investigation  has  served  to  establish  the  pattern  for  later  efforts  on  impfoved  theories 
for  study  of  free-edge  delamination  of  composite  laminates. 

The  equations  of  the  discrete  laminate  theory  are  well  known.  That  have  been 
included  in  this  report  for  completeness  and  for  the  purpose  of  pointing  out  the 
theoretical  assumptions  inherent  in  the  theory.  Mawenya  and  E)avis  [S]  had  previously 
presented  a  finite  element  implementation  of  the  theory,  but  they  provided  no  details 
and  did  not  apply  the  theory  to  free-edge  delamination. 

The  equations  of  the  discrete  laminate  theory,  essentially  treating  a  laminate  plate 
as  a  stacking  of  Mindlin  plates,  assumes  linear  variation  of  the  in-plane  displacements 
over  the  thickness  of  each  layer  ensuring  continuity  of  displacements  at  the  interfaces. 
As  part  of  the  present  research  effort,  the  equations  have  been  written  in  matrix  form 
such  that  the  matrix  of  operators  is  self-adjoint  in  an  appropriate  linear  vector  space. 
The  general  formulation  was  written  for  the  dynamic  problem  in  the  convolution 


122 


product  space  and  then,  for  the  elasiosiaiic  problem,  a  special  fonn  was  defined  in  the 
inner  product  space.  For  such  systems  of  equations,  standard  techniques  are  available 
for  the  construction  of  variational  principles  and  for  identification  of  consistent 

boundary  operators.  General  variational  formulations  with  extensions  and  useful 
1  specializations  have  been  explicitly  developed  for  the  problem.  One  specialized 

^  variational  formulation,  corresponding  to  the  popular  potential  energy  thcoryp  has  been 

used  to  develop  finite  element  procedures. 

Three  different  isoparametric  finite  element  interpolation  schemes  viz.,  the 
four-point  Lagrangian,  the  nine-point  [^grangian,  and  the  eight-point  serendipity,  have 
been  implemented  in  a  computer  program  written  initially  for  an  IBM  3081  mainframe 
computer  and  later  modified  to  run  on  a  CRAY-XMP/28.  These  finite  element 
procedures  were  verified  through  application  to  homogeneous  as  well  as  sandwich  plate 
problems  for  which  solutions  are  available.  Their  effectiveness  in  modelling  the 
stress-distribution  in  free-edge  delamination  specimens  has  b^n  examined. 

The  present  investigation  indicated  that  the  discrete  laminate  theory  of  laminated 
plates  is  quite  effective  in  modelling  displacements  in  plates  subjected  to  arbitrary 
transverse  loads.  The  shear  effects  can  be  allowed  for  satisfactorily.  For  such 
problems,  the  higher  order  elements  viz.,  the  nine-point  Lagrangian  and  the  eight-point 
serendipity,  performed  better  than  the  simple  four-point  Lagrangian  element  However, 
for  application  to  free-edge  deiamination,  the  entire  approach  is  inadequate.  The  stress 
distribution  obtained  for  the  example  problems  was  reasonably  gocxl  with  respect  to 
in-plane  stresses,  but  the  theory  could  not  give  reasonable  estimates  of  the  other  three 
components  of  stress.  The  traction-free  edge  condition  could  not  be  modelled.  The 
stresses  had  to  be  calculated  directly  from  the  stress-strain  relationships  because  use  of 

equilibrium  equations  for  determination  of  shear  stress  and  transverse  stress  would 
involve  numerical  differentiation  of  quantities  for  which  estimates  only  at  a  finite 
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number  of  points  were  available.  Mesh  refinement  to  get  a  sufficiently  large  number 
of  points  would  make  the  cost  of  analysis  prohibitive. 

Tlie  studies  showed  that  refinement  along  the  length  or  the  width  of  the  specimen 
had  relatively  little  effect  on  the  quality  of  results.  Refinement  over  the  thickness 
subdivision  of  each  lamina  into  a  number  of  sublayers.,  helped  improve  accuracy. 
This  suggests  that  a  distribution  of  in-plane  stresses  of  an  order  higher  than  linear 
over  the  thickness  of  each  lamina  might  represent  the  actual  stress  distribution  more 
closely. 

In  order  to  satisfy  the  traction-free  conditions  along  the  free-edges,  it  is  necessary 
that  edge  tractions  appear  as  field  variables  in  the  set  of  field  equations.  An 
alternative,  of  course,  is  to  use  Lagrange  multiplier  techniques  to  enforce  constraint* 

For  direct  use  of  equilibrium  equations  to  determine  the  shear  stresses  and  the 
direct  transverse  stress,  explicit  introduction  of  interface  tractions  as  field  variables  in 
the  theory  would  avoid  the  need  for  expensive  numerical  differentiations.  Also  this 
would  ensure  continuity  of  traction  across  interfaces  and  perhaps  yield  better 
approximation  for  the  interfacial  stresses. 

The  discrete  laminate  theory  discussed  in  this  report  is  based  on  assumptions 
regarding  transverse  and  in-plane  displacements.  An  alternative  is  to  assume  variation 
of  in-plane  stresses  and  to  derive  the  other  stress  components  through  equilibrium 
equations.  If  force  resultants  appear  in  the  expressions  for  stresses,  and  are  regarded  as 
field  variables  of  the  problem,  the  constitutive  relationships  for  these  need  to  be 
established.  If  it  is  assumed  that  there  is  no  interfacial  slip,  it  would  be  impossible 
for  any  layer  to  deform  independently  of  the  others.  This  would  necessarily  lead  to  a 
coupling  in  the  constitutive  relations  for  the  force  resultants  of  individual  layers. 

Pagano’s  [lO]  theory  which  uses  the  assumption  of  linear  variation  of  in-plane 
stress  components  over  the  thickness  of  each  layer  or  sublayer,  satisfies  equilibrium 
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pointwise,  satisfies  constitutive  equations  and  inierfactal  continuity  of  tractions  as  well 
as  displacements,  and  can  satisfy  traction  boundary  conditions  for  free-edge  delamination 
specimens  exactly,  appears  to  be  an  appropriate  approach  for  determination  of  stress 
fields  in  composite  laminate  plates  with  natural  boundary  conditions.  The  case  of 
free-edge  delamination  specimens  is  a  specialization  of  the  general  theory.  The  theory 
has  been  difficult  to  use  because  of  the  large  number  of  field  variables  involved  and 
limitations  on  computational  capabilities.  A  research  effort  directed  towards 
development  of  finite  element  models  based  on  Pagano's  theory  or  development  of  other 
niethcxls  of  solution  of  the  .set  of  differential  equations  could  be  useful. 


j 

a 


125 


REFERENCES 


t 

1.  Srinivas,  S.,  A  Refined  Analysis  of  Composite  Laminates,  J.  Sound  &  Vib^ 

Vol.  30.  No.  4,  495-507.  1973. 

2.  Sun.  C.T.  and  Whitney.  J.M.,  Theories  for  the  Dynamic  Response  of  Laminat¬ 
ed  Plates,  AIAA  J..  Vol.  11.  178-183,  1973. 

3.  Reddy.  J.  N.,  A  Penalty  Plate-Mending  Plement  for  the  Analysis  of  Laminat¬ 
ed  Anisotropic  Co?nposite  Plates,  Int.  J.  Numer.  Methods  in  Engrg..  Vol.  15,  pp. 

1187-1206,  1980. 

4.  Reddy,  J,  Dynamic  (Transient)  Analysis  of  Layered  Anisotropic  Composite- 
Material  Plates,  Int.  J.  Numer.  Methods  in  Engrg.,  Vol.  19,  237-255,  1983. 

5.  Mawenya,  A.S.  and  Davis,  JD.,  Finite  element  bending  analysis  of  Finite  Ele¬ 
ment  Bending  Analysis  of  Multilayer  Plates,  Int  J.  Numer.  Methods  in  Engrg.. 

Vol.  8,  215-225.  1974. 

6.  Sandhu,  R.  S.,  and  Pister,  K.  S.,  A  Variational  Principle  for  Linear  Coupled 
Field  Problems,  Int  J.  Eng.  Sci.,  Vol.  8,  986-999,  1970. 

7.  Sandhu,  R.  S.,  and  Pister,  K.  S.,  Variational  Princip}les  for  Boundary  Value 
and  Initial  Boundary  Value  Problems,  Int  J.  Solids  Struct,  Vol.  7,  639-654, 

1971. 

8.  Sandhu,  R.  S.,  and  Salaam,  D.,  Variational  Formulation  of  Linear  Problems  ) 

vHth  Nonhomogenous  Boundary  Conditions  and  Internal  Discontinuities.  Comp. 

J 

Meth.  Appl.  Mech.  Eng^  Vol.  7,  75-91,  1975. 

9.  Sandhu,  R.  S.,  Variational  Principles  for  Finite  Element  Approximations  in 
Finite  Elements  in  Water  Resources  Engineering,  G.  F.  Pinder  and  C.  A.  Breb- 
bia  (eds.),  Peniech  Press.  1976. 


126 


10.  Pagano  NJ.,  Stress  Fields  in  Composite  Ixtminaies,  Int.  J.  Solids  Structures, 
Voi.  14,  385-400,  1978. 

11.  Chang,  C.  C.,  Finite  Element  Analysis  of  Laminated  Composite  Free-Edge 
Delanunation  Specimens,  Ph.D.  Dissertation,  TTie  Ohio  State  University,  1987. 

12.  Sandhu,  R.  S^  Sierakowski,  R.  U  Wolfe,  W.  E,  and  Dandan,  R.  A,  Finite  Ele¬ 
ment  Analysis  of  Laminated  Composite  Axisymmetric  Solids.  Vol.  I  -  Theory 
and  Applications,  Wright  Laboratory  Technical  Report  WL-TR-9 1-3020,  Vol.  1, 
Wrighi-Patterson  Air  Force  Base.  Ohio,  1991. 

13.  Mikhlin,  S.  C.,  7'he  Problem  of  the  Minimum  of  a  Quadratic  Functional, 
ilolden-Day,  San  l-'rancisco,  1965. 

14.  Plantema,  FJ.,  Sandwich  Construction,  John  Wiley  and  Sons,  1966. 

15.  Khatua  T.P  and  Cheung  Y.K,  Bending  and  Vibration  of  MtdtUayer  Sandwich 
Beams  and  Plates,  Inl.  J.  Numer.  Methods  in  Engrg.,  Vol.  6,  11-24,  1973. 

16.  Pagano  NJ.  and  Hatfield  SJ,  Elastic  Behaviour  of  Multilayered  Bidirectional 
Sandwich  Beam  and  Plates,  AlAA  J,  Vol.  10,  931-933,  1972. 

17.  Pagano  NJ.,  Exact  Solution  for  Rectangular  Bidirectional  Composites  and 
Sandwich  Plates,  J.  Composite  Matls.,  Vol.  4,  20-24,  1970. 

18.  Zienkiewicz,  O.C.,  The  Finite  Element  Method,  3rd  edition,  McGraw-Hill,  1977. 

19.  Ugural,  A.Ch  Stresses  in  Plates  and  Shells,  McGraw-Hill,  1981. 

20.  .^zar,  J.  J.,  Bending  theory  of  multilayer  sandwich  pltdes,  AlAA  J.,  Vol.  6, 
2166-2169,  1969. 


127 


Appendix  A 

VARIATIONAL  FORMULATION 


Often,  obtaining  an  approximate  solution  to  a  coupled  boundary  value  problem 
relies  on  appropriate  variational  formulation.  Following  Sandhu's  (6]  (?]  [8]  [9], 
extension  of  Mikhlin's  [13]  basic  variational  theorem  to  coupled  linear  boundary  value 
problems  including  nonhomogenous  boundary  condition,  we  present  here  a  summary  of 
the  basic  concepts  for  setting  up  the  variational  formulation  applicable  to  the  problem 
of  laminated  plates. 

A.1  PRELIMINARIES 

A.1.1  Boundary  Value  Problem 

Consider  the  boundary  value  problem 


(A.l) 


Au  =  /  on  R 


(A.2) 


Cu  ~  g  on  QR 


where  QJ?  is  the  boundary  of  the  open  connected  region  R  in  an  euclidean  space.  S  is 
the  closure  of  R.  A  and  C  are  linear  bounded  operators.  Let  Vj,  and  be  linear 
vector  spaces  defined  on  the  regions  indicated  by  the  subscripts,  and  be  dense 

subsets  in  V;,  and  respectively.  Then  the  differential  operators  A  and  C  can  be 
regarded  as  the  transformations 


CA.3) 
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A.1.2  Bilinear  Mapping 

Let  V  and  S  be  linear  vector  spaces.  A  bilinear  mapping  fi:VjcV-»S  assigns  to 
each  ordered  pair  of  vectors  w.vgV  an  element  in  S.  Furthermore,  bilinearity  is 
satisfied  for  u,  ,Uj,  v, ,  v^.u ,  v  €  V,  if 


Kau,  +«2,  v)  =  aB(u^ ,  v)  +  B(u^,v) 

(A.4) 

4 

SCu.oVj  +  Vj)  =  aB(u,v^)  + 

where  a  is  scalar.  For  convenience,  we  shall  use  Die  notation. 

(A,5) 

£jj(«,v)  =  <  u.v 

CA.6) 

To  set  up  a  variationiii  formulation,  symmetric,  nondegenerate  bilinear 

used,  t.eH 

mappings  are 

<U,V>  JJ= 

and 

(A.7) 

<u,v>  =  0  for  all  V  if  and  only  if  u  =  () 

A.  1.3  Self 'Ad joint  Operator 

(A.8) 

An  operator  A*  on  V  is  said  to  be  the  adjoint  of  A  with  respect 

bilinear  mapping  Bg'.VxV S.  where  S  is  a  linear  vector  space,  if 

to  symmetric 

<«.Av>^=  <v.A*u>j,  +  £i^(v,h) 

(A.9) 

for  all  u  and  v  6  V  and  where  Z)j,^(M,v)  represents  quantities  associated  with 

boundary  ^  of  R.  If  A-  A*,  then  A  is  said  to  be  self-adjoint.  If  A  is 

operator,  then  Z)gj,(v,u)  is  antisymmetric,  ie., 

a  self-adjoint 

J 

D^(v.u)  = 

CA.IO) 

Furthermore,  A  is  said  to  be  symmetric  with  respect  to  the  bilinear  mapping,  if 

<  «,Av>j^  =  <  v,  Au  >jp  (A.ll) 

The  boundary  operator  ('  is  said  to  be  consistent  with  the  self-adjoint  operator  A  if 
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D^/v,u)=  <u,Cv>.^^~<v,Cu>^^ 


(A.12) 


A.1.4  Gateaux  Differentia) 

If  fl;V“*5,  where  V  is  siati  i)iat  if  u,  Q€V',  ih  Aii€V’  for  scalar  X, 
the  Gateaux  differential  of  fl  (u)  along  a  path  ii  is  defined  by 

S.  n(u)  =  lim  +  (A.] 3) 

“  x-o  X 

where  u  is  referred  in  as  tlie  patli. 

A.2  THEOREM 

I'or  the  field  equations  (A.l)  we  del'ine 

n(«)  =  <w ,  A«>  j(,  ~2  <u  ,/>  jj  (A.H) 

The  Gateaux  differential  of  fl  is: 

8  Sliu)  -lim  ^  +  <i/4-XG,/>-<u,Au>  +  2  <u,f>  (a. 15) 

"  x-o  X 

=  <  u,Au>  +  <u,Au>— 2  <«,/> 

=  2  <  u ,  Au  —  f> 

The  Gateaux  differential  vanishes  at  the  solution  u  =  u,  where  Au,  — /  =  0.  Conversly, 
if  SgfKu)  vanishes  for  all  u,  nondegeneracy  of  <  ,  >  implies  Au,  — /  =  0  If  the  range 
of  the  bilinear  mapping  is  the  real  line,  vanishing  of  the  function  fl  would  imply  its 
minimum,  maximum,  or  stationary  value,  depending  upon  the  operator  A  being  positive, 
negative  or  semi-definite. 


f 

V 


J 


13(1 


A.3  UNEAR  COUPLED  PROBLEMS 


The  above  discussion  for  a  si nj>le- valued  function  u  can  be  extended  to  the  case  of 
several  variables.  If  there  are  n  variables,  V  is  defined  as  the  direct  sum 

V  =  V,  +  V2  + - +  (A.16) 

and  an  element  uCV  is  an  n-tuple  (w,  ,u, . ,«„)  with  «,  €  for  i-1,  2,  — ,  n.  A 

bilinear  mapping  on  V  is  defined  as 

<u,v>  =  <u,,v,>  ,+  <u,,V2>  . +  <w^^.v^>  CA.17) 

where  <  .  >^  is  defined  for  components  u, ,v,  of  {u,},{v,}  respectively. 

If  the  I'ield  anti  boundary  condition  of  a  linear  coupled  boundary  value  problem 
are; 


n 

^C.jUj  =  g.  on  Q/t  i=1,2,...^  (A.19) 

j-i 

the  governing  functional  based  on  Fqs.  (A.  18)  and  (A.  19)  is 

n  n  if  n 

n  Cu)  =  £  <  u, .  £  A  -  2/,  > ;,  +  £  <  «, .  £  C.  «  ,  -  2g.  >  ^  (A.20) 

i-i  J-I  i-i  y-1 

The  set  of  operators  A^j  is  said  lo  be  self-adjoint  with  respect  to  the  bilinear  mapping, 
if 


ft 


<  V 


T  A  M.>  ^  =  <  V  a  ,A^' 

'  ^  u  j  ^  /  St 


(AJl) 


yi 

Where  represents  quantities  associated  with  boundary  QR  of  R.  The  boundary 

operators  are  said  to  be  consistent  with  the  field  operator  A,j  if 


D  lu  ,  V.)  =  V  <  ti  ,C  V,  >  „  —  <  V  ,  T*  C, ,«  ,  > 

^P  j  t  ^  j  yi  I  0^  tj  j 


BP 


(A.22) 


Appendix  B 

SOLUTION  OF  SANDWICH  PLATE 


B.1  PRELIMINARIES 

liiised  on  Pliintemit  iM],  ihe  series  solmion  Tor  sandwich  plate  was  calculated.  In 
this  case,  tlie  Ijcndin^  ol  plates  is  assumed  to  '>e  due  to  bending  of  stiff  layers  and  the 
shear  deformation  of  the  core  layer,  so  that 


w  =  +  w 

h  i 


tB.l) 


where  the  transverse  load  is  the  only  applied  load  and  the  stiff  layers  are  isotropic, 
satisfies  the  equation: 


Z)  V  w  =  9 


(B.2) 


where 


D  = 


Et\ 


12(1 -v) 

and  w,  satisfies  the  equation: 


-5  V  w,  =  q 


(B.3) 


where 


S  = 


 2 


t. 


— G, 


i 


1.12 


B.2  METHOD  OF  SOLUTION 

For  a  simply  supported  rectangular  plate  the  transverse  displacement  due  to 
bending  may  be  represented  by  Fourier  series  in  the  form: 

w(x,y)  =  £  a^sin  sin-Sp^  CB.4) 

m*-  !  n-l 

and  corresponding  load  by: 

=££«„, sin  ip  (B^) 

HI- 1  n- 1  °  ^ 

*  *  * 

Multiplying  both  sides  of  (H-5)  by  sin  ^  sin  ^  and  integrating  over  the  domain, 

Q  b 

for  q(x,y)=  q,^ 


16a 

a  =  ^ 

2 

n  mn 

Substituting  (B.4)  and  (B^)  into  (BJ)  and  evaluating  v/Jix,y)  is: 


££«in(mw/risin(n.ry/»  , ,3,, 

“IT  D  „  „  mn[(m/a)  +  Cn/6)^]^ 

Similarly,  (B.3)  can  be  solved  to  yield: 


Ca6) 


,  ^  y*  sin  (inTrx/o)sin  (ntry/fe) 

_  „  m»l(.m/a)^  +in/b)^] 


The  series  solution  was  obtained  through  summation  of  151  terms  in  the  case  of 
isotropic  sandwich  plate  for  a-l^lO  inches  and  q^=  i.O  Ih/in^,  In  the  orthotropic  case 
the  results  given  by  Mawenya  and  Davis  [5]  were  used  for  comparison  purposes. 
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